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ABSTRACT 



We extend a systematic renormalization procedure for quantum field theory to in- 
clude particle masses and present several applications. We use a Hamiltonian formu- 
lation and light-front quantization because this may produce a convergent Fock-space 
expansion. The QCD Hamiltonian is systematically renormalized to second order in 
the strong coupling and the Fock-space expansion is truncated to lowest order to 
produce a finite-dimensional Hamiltonian matrix. The renormalized Hamiltonian is 
used to calculate the spectra of the bb and cc mesons as a lowest-order test of our 
procedure for full QCD. 

The analytic determination of the renormalized Hamiltonian matrix generates 
expressions that must be numerically integrated to generate quantitative results. The 
efficiency of the numerical calculation depends on how well the basis functions can 
approximate the real state. We examine the effectiveness of using Basis-Splines (B- 
Sphnes) to represent QCD states. After briefly describing these functions, we test 
them using the one- and two-dimensional harmonic oscillator problems. We test their 
ability to represent realistic wavefunctions by using them to find the glueball mass 
spectrum. 

An efficient algorithm for numerically calculating the matrix elements in the glue- 
ball and meson problems is necessary because the calculation is numerically intensive. 
We describe our algorithm and discuss its parallel-cpu implementation. 
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CHAPTER 1 



Introduction 



The physics of quarks and gluons is described by quantum chromodynamics (QCD), 
which is universally believed to be the fundamental theory of the strong interaction. 
High energy QCD systems (small spatial separations) behave as though the particles 
interact weakly (asymptotic freedom). Low energy systems (large spatial separations) 
interact so strongly that it is not possible to isolate a quark or gluon (confinement). 
The physics in the confining and asymptotically free systems is so different that it 
is has not yet been possible to write down a solution to QCD that is valid in both 
regimes. Historically, difficult theories like QCD are first solved for simple systems. 
Solving these systems provides an initial test of the theory before more complex sys- 
tems are attacked and a development ground for new techniques. 

Heavy quark bound states are relatively simple because the large quark mass 
dominates the spectrum and the average separation of quarks and gluons is small, 
so that dynamics are largely described by asymptotic freedom. This also allows the 
interactions to be treated nonrelativistically |, |, ||. Techniques that have 
been applied to the heavy quark system include nonrelativistic potential models and 
nonrelativistic QCD (NRQCD). Potential models p, obtain good predictive power 
by approximating the interaction between particles using nonrelativistic potentials, 
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although these potentials are not directly obtained from QCD. Predictions include 
the T/T' splitting ||^, spin-splitting, transition rates and hadronic transitions 0. 

NRQCD uses a nonrelativistic expansion of QCD operators to generate an approx- 
imate Lagrangian that can be systematically improved 0. Heavy-quark spectroscopy 
can be determined on the lattice |TD|, which can then be used to determine the strong 
coupling constant, which is important because the coupling can be compared to per- 
turbative results^ and could reveal the existence of new physics fl^ . Both potential 
models and NRQCD have proven to be predictive theories for nonrelativistic systems. 
Despite both theories being improvable, a relativistic treatment of QCD is necessary 
to treat systems of lighter quarks. 

There are many problems one must solve (or circumvent) to solve full QCD. In 
Euclidean or equal-time coordinates vacuum fluctuations prevent using a constituent 
picture for low energy hadrons. Also, since QCD is a field theory, there are divergences 
that must be regulated, often breaking symmetries. Any solution to QCD will require 
the use of renormalization to remove the dependence on the regulator and to restore 
broken symmetries. 

Hadrons are composed of at least two (meson) or three (baryon) quarks.0 The 
quarks are bound by the strong force through the exchange of color-charged gluons, 
which are in turn bound together by the same mechanism. Hadronic bound states 
(and glueball states) are color singlets; there is no net color. As quarks in a baryon 
are separated, the energy of the system increases, eventually allowing the creation 

of new quark-antiquark pairs. These pairs form bound states with the separated 

^There are also perturbative models of NRQCD. See pH |. 

^In this overview we do not distinguish between quarks and antiquaries. 
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quarks, preventing the isolation of any individual quark. If we consider low-energy 
bound states, the particles interact weakly and we may be able to make a constituent 
approximation. However in a standard equal-time picture Heisenberg's uncertainty 
principle allows quark-antiquark pairs to fluctuate out of the vacuum. These quarks 
can then interact with any hadronic bound states, preventing the state from being 
successfully represented by two or three quarks. Moreover, these vacuum fluctuations 
are usually thought to be critical for strong interactions. 

We would like to derive a constituent picture for mesons from QCD, with mesonic 
states and masses determined by a Schrodinger equation: 

H\^) = E\^), (1.1) 

where we can approximate the eigenstate with a truncated Fock-space expansion, 

1^) = <Pgg\(l(l) + 4>q<lg\qq9) + " " •• (1-2) 

E is the energy eigenvalue of the state and we use shorthand notation for the Fock 
states where g is a quark, q an antiquark, and (7 is a gluon. Unfortunately a Fock-space 
expansion is not usually thought to be practical, because the complicated equal-time 
vacuum forces bound-states to contain an infinite number of particles that are part 
of the physical vacuum on which hadrons are built. It may be possible, however, 
to obtain a convergent Fock-space expansion if we work in light-front coordinates 
because the free energy of a state increases at least like the number of particles 



squared (Sec. p.l|) . This means all states with many particles are high-energy states. 



Inspired by the work of Dyson [|l^ , Wilson [|T^ , Glazek and Wilson [T^ , and Wegner 



16| , significant work has been done to perturbatively derive light-front Hamiltonians 
in the full Fock-space, neglecting zero modes |, |, |, [l^, |18|, |l|, |g, 0, |2|, |3|, |4 . 
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Local interactions in a field theory couple states with arbitrarily large differences 
in free energy, which will invalidate our Fock-space expansion if we can not decouple 
high-energy states from low-energy states. One way we can force high-energy states to 
decouple from low-energy states is to simply truncate the Fock-space expansion. The 
truncated states can be important, and their removal may discard important physics. 
We can retain all of the physics if we place a cutoff on change in invariant mass 
between states, and then use a similarity transformation ||15|, [16| to map the effects 
from mixing of low- and high-energy states to interactions involving only low-energy 
states. 

The importance of high-energy states can be appreciated by considering the second- 
order perturbation theory correction to a free energy: 



N 





l< 


m 


V 




{m 


Ho 


m) - 


- {n 


Ho 


n) 



= E /^i^i ! /iWri^v - (1-3) 



The sum can diverge as — > oo if the matrix elements between states with large 
differences in energy do not fall off sufficiently rapidly. We regularize the theory 
by implementing a gaussian cutoff on interactions between states with differing free 
invariant masses^. This cutoff serves two purposes. It makes perturbative corrections 
like Eq. ( p.. 3D finite and it reduces the coupling between states with large differences 
in invariant mass, which further helps justify our Fock-space expansion by reducing 
the importance of high-energy physics for low-energy states. 

The cutoff we use violates Lorentz covariance and gauge invariance, so we can 

not renormalize the Hamiltonian exclusively through the redefinition of masses and 

•^The invariant mass can be thought of as the energy of the state after subtracting the energy 
due to translations. 
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canonical couplings. Renormalization must be completed by requiring the Hamilto- 
nian to produce cutoff-independent physical quantities and by requiring it to obey 
the physical principles of the theory that are not violated by the cutoff. These re- 
quirements, and the assumption that this can be done perturbatively, are sufficient to 
determine the Hamiltonian so that it will give results consistent with all the physical 
principles of the theory, even those violated by the cutoff. The most powerful char- 
acteristic of this approach is that it systematically "repairs" the theory and requires 
only the fundamental parameters of the canonical Hamiltonian. 



Recently, Allen [^S], ^ derived recursion relations that systematically deter- 
mine the Hamiltonian for massless theories order by order in the coupling and he 
applied this method to pure-glue QCD. In this dissertation we extend the recursion 
relations to full QCD and present three applications. We first show how matrix 
elements of the renormalized Hamiltonian are calculated with our extension of the 
method, which has been published in |^ . Secondly we apply new numerical methods 
to verify Allen's previous glueball spectrum calculation |^B|, . The final application 
uses these numerical methods to calculate heavy quarkonia {qq) spectra. 

In light-front field theory the Hamiltonian is trivially related to the invariant 
mass operator (IMO). The IMO is a natural operator to use because it is boost and 
rotationally invariant (although rotational invariance is realized dynamically). We 
calculate the IMO in this dissertation, but refer to it as the Hamiltonian throughout 
since it governs the time evolution of all states. 

The Hamiltonian matrix is calculated analytically to 0{as) in a plane- wave basis. 
We then calculate the matrix elements in a basis we choose to approximate the real 
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meson states (i.e., eigenstates of the Hamiltonian) . The real meson states are approxi- 
mated by the first term in a Fock-space expansion, quark-antiquark pairs. Solving the 
eigensystem gives the invariant masses (eigenvalues) of the system and the coefficients 
(eigenvectors) needed to build the eigenstates from the basis functions. 

The numerical calculation of the Hamiltonian matrix requires us to choose a set of 
basis functions to approximate the real states. We use Basis-splines (B-splines) to rep- 
resent the longitudinal and transverse momentum states. Although their derivation 
from mathematical first principles is non-trivial, their benefit is easily understood. 
First, B-splines do not oscillate, so they do not add extra oscillations to the function 
that must be integrated to determine each matrix element .0 Second, each spline has 
only a limited spatial overlap with the other splines; thus, there can be a significant 
number of matrix elements that vanish if the interactions are local because the basis 
functions do not overlap. 

The most substantial drawback to using B-splines is that they are non-orthogonal. 
This means that an overlap matrix must be calculated to solve Schrodinger's equation. 
The main inconvenience is that a general eigensystem problem must be solved and 
efficient methods that compute only low-lying eigenvalues can not be used. 

The numerical calculation of the Hamiltonian matrix is very cpu-intensive because 
each matrix element includes a five- dimensional integral. Unfortunately one can not 
know ahead of time how accurately each matrix element must be calculated. For 
instance, if all matrix elements are calculated to the same precision, much computer 

time is wasted because all elements are not equally important when determining 

^The matrix elements are approximated using Monte-Carlo integration. Monte-Carlo is a pow- 
erful method to approximate many-dimensional integrals, but it does not work well with oscillating 
integrands. 
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the low-lying eigenvalues. With no unambiguous method to determine the relative 
importance of various matrix elements, there is also no obvious algorithm to solve the 
problem efficiently. Fortunately each matrix element can be calculated independently, 
so most algorithms can be executed in a parallel computing environment. 

In Chapter ^ we discuss our motivation for using light-front coordinates and our 
choice of regularization and renormalization techniques. We then review a recent 
technique that builds on previous work, systematically determining the renormalized 
Hamiltonian for massless theories, and discuss how this method is generalized to 
theories with massive particles. We conclude Chapter ^ by listing which terms in 
the Hamiltonian are needed for an 0{a) calculation, followed by a discussion of how 
these terms are combined to explicitly cancel light-front divergences. 

We present B-splines and show some simple applications in Chapter 0. In Chapter 
H we introduce the basis used to represent the real meson states. This basis is tested 
by computing the glueball mass spectrum in Section . 

The matrix elements of the Hamiltonian are presented in Chapters ^, |^, and |^ 
Chapter ^ presents some matrix elements of the full Hamiltonian in a plane wave 
basis. The matrix elements of the second-order renormalized Hamiltonian in a plane 
wave basis are derived in Chapter ^. Finally the matrix elements of the Hamiltonian 
in the approximate-state basis are given in Chapter |^. 

Although there is a lot of analytic work required to solve this problem, quantitative 
results can only be obtained after applying numerical methods to attack the problem. 
We discuss some of the numerical considerations when deriving the matrix elements 
and describe the general algorithm and how it is parallelized in Chapter 0. 
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Finally, we present our analysis of our numerical methods and the calculation of 
the meson spectrum in Chapter |^. We make concluding remarks in Chapter 

This dissertation builds upon previous work and applies it to a familiar problem 
using numerical techniques that may be similar to other calculations unknown to the 
author. We believe the following parts of this calculation are original work: 

• We included masses in Renormalization Group calculations using the method 



in 18 . 



Calculation of the renormalized light-front QCD Hamiltonian to 0{a) using the 



renormalization methods from ||T^, |T8|, ^| 



Use of B-spline basis function to calculate the glueball mass spectrum. 

Numerical calculation of the heavy meson mass spectrum for our renormalized 
QCD Hamiltonian. 

Development of an efficient algorithm to calculate the Hamiltonian matrix. 



8 



CHAPTER 2 



Motivations and Formalism Overview 



In this chapter we motivate and give an overview of our method. This includes 
a discussion of hght-front coordinates as well as discussions about regularization and 
renormalization. We then discuss a recent development that systematically deter- 
mines the renormalized Hamiltonian by requiring it to produce cutoff-independent 
results and to obey all physical principles not violated by the cutoff |^ |2^, P7 |. 



We then review the generalization to include massive particles from p8|. In Section 
|2.6| we list all of the diagrams that contribute to the second-order Hamiltonian. Fi- 
nally we describe how these diagrams are combined to explicitly cancel the light-front 
divergences. 

This chapter is intended to give a good understanding of the entire approach 
without saturating the reader with too many details, which can be found in the 
references. 

2.1 Light-Front Coordinates 

Light-front coordinates may be considered a rotation of equal-time coordinates 



such that light-front time is along the forward light-cone (See Fig |2.1| ). Explicitly we 
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X 








Figure 2.1: Light-from coordinates in 3 + 1 dimensions. x+ and the light- front 
time and longitudinal space components are a linear combination of the equal-time 
3-direction and time. 



write the light-front time coordinate: 

a;+ = x° + x^ (2.1) 
and the light-front longitudinal space coordinate as: 

x-=x^-x^. (2.2) 

The scalar product of light-front vectors is: 

a ■ b = -a~^b~ + -a~b'^ — a± ■ h±, (2.3) 

where and are the transverse components. This shows that the light-front 
energy, which is conjugate to light-front time, is p~ and the longitudinal momentum 
is p^. An on- mass-shell particle obeys > p^ which implies p^ > 0. Finally a free 
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particle of mass m has the hght-front energy 

r = (2.4) 

The hght-front energy is finite unless the longitudinal momentum is zero. The study 
of these infinite-energy zero modes is rich.0 Since they are infinite energy modes we 
discard them in our approach, although their effects will need to be added to the 
Hamiltonian. Using the light-front dispersion relation we can write the free energy of 

Qb sIjciIjG clS! 

1 N 2 I 2 

^_^^^p^_ (2.5) 

where is the total longitudinal momentum of the state, N is the number of particles 
in the state, the index i refers to each particle in the state, and Xj is the i^^ particle's 
longitudinal momentum fraction. If we assume all particles in the state have similar 
masses and transverse momenta, we can write the energy as: 

p_ _ (pj + m^)typicai ^ 1 ^ ^2.6) 

The free energy is minimized if the particles equally share the longitudinal momentum 
{xi = This means the minimum free energy of a light-front state behaves like: 

^min. ~Ar2(p^ + m2),ypi,,i. (2.7) 

Since the energy of the state increases like the number of particles squared, a Fock- 
space expansion may be justified. 

See appendix ^ for a complete presentation of the light-front formalism used in 
this paper. 

^Reference is a review of theories on the hght-front that includes a discussion of zero-modes 
and many references to zero-mode calculations. 
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2.2 Regular izat ion 

In the next two sections we discuss the procedure we use to regulate and renor- 



mahze the Hamiltonian. Our discussion follows the one in [0 and we do not cite 
each instance an argument from this reference is used. 

If particle number is not conserved in a theory, there can be arbitrarily many 
particles in a state. In addition each state can be coupled to every other state through 
a matrix element of the Hamiltonian, yielding a Hamiltonian with an infinite number 
of matrix elements. We must regularize the Hamiltonian so it is finite dimensional. 
We could simply limit the number of particles in a state. This type of truncation is 



called a Tamm-Dancoff truncation |3T], |3^ and suffers from divergent sensitivity 



to the precise form of the truncation, and removing this sensitivity can be as hard as 
solving the full theory. 

Consider a Hamiltonian that can be written as; 

H = Ho + V, (2.8) 

where V can be considered a perturbation. To second order in perturbation theory, 
the energy of a state |n), is given by: 

E„ = {n \Ho\ n) + {n \V\ n) - ^ , .JM^'fi^i ^ - (2-9) 

n^n {mlHolm) - {n\Ho\n} 

If this sum is over a finite number of states, and the matrix elements are finite, the 
sum will be finite. However, if the sum is over an infinite number of states, the sum 
may not converge if the matrix elements between states do not fall off fast enough. 
One way we can ensure this sum is finite is to place a cutoff on the Hamiltonian that 
reduces the coupling between states with arbitrarily large differences in free energy. 
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Since perturbative calculations can only lead to reasonable results when the cutoff 
is a function of the free energy variables and we want a cutoff that will also prevent 
small energy denominators, we choose to place a smooth cutoff on change in invariant- 
mass (free energy) between the final and initial states [Eq. ( ^.17] )]. Like all cutoffs, 
this violates Lorentz covariance and gauge invariance. However, it is invariant under 
boosts and rotations about the 3-axis. 

The light-front dispersion relation shows that the energy of a state rapidly in- 
creases with the number of particles, and the cutoff reduces the coupling between 
high- and low-energy states. These two facts support the use of a Fock-space expan- 
sion. 

2.3 Renormalization 

We regulate the Hamiltonian by applying a smooth cutoff on change in invariant- 
mass between the final and initial states. This introduces cutoff dependence in the 
Hamiltonian, and it breaks Lorentz covariance and gauge invariance. We restrict the 
form of the Hamiltonian by requiring it to obey all physical principles unviolated by 
the cutoff and use the similarity renormalization scheme to perturbatively remove 
the cutoff dependence. An 0{a) perturbative renormalization of the regularized 
Hamiltonian does three things. It removes cutoff dependence and restores broken 
symmetries to 0{a). 

For notational consistency]^ with Ref. (which we follow in this section) we 
remind the reader that if we neglect the total transverse momentum of a system,Q the 



6 



Some of the notation is this section is shorthand to simphfy the discussion. The fuU notation is 



introduced in Section 2.4 



^Transverse boosts are kinematic on the hght-front, so aU of the interactions in the Hamiltonian 
are in the IMO. 
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Hamiltonian is trivially related to the Invariant Mass Operator (IMO) by: 

_ (2.10) 
We split the IMO into free and interacting parts: 

M\A) = Ml,, + MUA)- (2.11) 

The IMO will produce cutoff independent results if it is unitarily equivalent to itself 
with an infinite cutoff and satisfies the relation:^ 

M^{A) = U{A, A')M\A')U^{A, A'), (2.12) 

where f/ is a unitary operator that reduces the cutoff from A' to A. The transforma- 



tion, which is a simplified version of the one introduced by Wegner |T6[, is determined 
by the differential equation: 

dU{A,A') 



d(A-4) 

with the boundary condition: 



r(A)f/(A,A'), (2.13) 



[/(A,A) = 1. (2.14) 

We treat the free part of the IMO (A^lee) nonperturbatively and use: 

T{A)=[Ml,,,M'{A)]. (2.15) 

If we can determine T(A), and therefore U{A,A'), then we can use the perturbative 
transformation repeatedly to lower the cutoff as long as the couplings are not too large. 

We can choose the initial cutoff to be as large as we want assuming the couplings 

®We assume the A oo limit exists. 
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in the large cutoff limit do not grow, which is true in an asymptotically free theory. 
This implies that M.'^{h) is unitarily equivalent to limA^oo -^^(A), the IMO with no 
cutoff. Thus 7W^(A) will give cutoff-independent results. 

We continue with the discussion of the exact form of the transformation and 
how we add restrictions from unviolated physical principles to fix the form of the 



renormalized Hamiltonian after more notation has been introduced in Section |2^ . 
2.4 Review of the Systematic Approach 



In this section we introduce more of the notation developed in Refs. ||2^, ^ ^| 
and outline the method to systematically determine the renormalized Hamiltonian. 
Formalism that is necessary for a detailed understanding of this method but that we 
do not repeat in this paper can be found in this earlier work. 

We want to find the regulated invariant-mass operator, A1^((yf^, m. A), which is 
trivially related to the Hamiltonian. It can be split into a free part (which contains 
implicit mass dependence) and an interacting part: 

M\g^,m,k) = Ml,,{m) + Ml^,{g^,m,K). (2.16) 

Since the method treats J^f^^{gj^,m, A) perturbatively, we put the particle-mass term 
in A1|gg(m), to treat it non-perturbatively; however, M.f^^[gj^^m, K) will still have 
mass dependence. The matrix elements of M.'^{gj^,m, h) are written 

{F\M\g,,m,k)\I) = {F\MlM\I) + {F\MU9.,^,m) 

= MUF\I) + e-^{F\V{g„m,A)\I), (2.17) 
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where \F) and |/) are eigenstates of the free invariant-mass operator with eigenvalues 
Mp and M|, and Apj is the difference of these eigenvalues. V{g^,m, A) is the inter- 
acting part of the invariant-mass operator with the Gaussian cutoff factor removed 
and is called the "reduced interaction." The Gaussian cutoff on change in invariant- 
mass is what regulates the theory (Sec. p^.2| ) and it has no effect on the free part of 
the Hamiltonian. 

We expand V{g^, m, A) in powers of the running coupling, g^: 

oo 

V{g,,m,A) = Y.glV^^\m,A), (2.18) 

r=l 

where V^^"^ is the canonical interaction and the y^^-^^m, A)'s are non-canonical in- 
teractions in the scalar theory. These non-canonical operators can be thought of as 
counterterms in a traditional approach. Note that (/^ implicitly depends on m. The 
assumption that this expansion exists is equivalent to a restricted coupling coher- 



ence . We now continue the discussion of the unitary transformation described in 
Sec. 



Allen shows in that the perturbative version of the unitary transformation 



that lowers the cutoff [derived from Eq. (|2.12|) ] is: 



V{A) - V{A') = SV, (2.19) 

where SV is the change in the reduced reaction and is a function of both A and A'. 
Since the free part of the Hamiltonian is independent of the cutoff, calculating 6V 
gives us the change in the Hamiltonian when lowering the cutoff from A' to A. The 
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matrix elements of SV are found to beQ 

{F\ 5V \I) =W {F\ n^') \K) (K\ n^') \I) Ti'^'^'XF, K, I) 

^ K 

^\ E (^1 ^(A') 1^) (^1 ^(A') 1^) (^1 ^(A') |/) Tf ^'^(F, L, /) 

./m4 



+0 ([l^(A')]'j , (2.20) 
where the sums are over complete sets of states. Tg^'^ ■* (F, I) and T^^'^ ■* (F, K, L, I) 



are: 



Tf'^')(F,ir,/) = - (e^^'-^^-xA.. - e^^"'^--^-^) , (2.21) 

\^FK ^KlJ 



and 



Tf'^')(F,i^,L,/): 
/ 1 1 



.^KL ^LI / \^KI ^FK , 

2A'— lA A 



AkL ^LiJ K^KL^LI + ^FK^KlJ 

KL^LI> 

' 1 1 \ / 1 1 \ 



A;.^ Ar,j, \Ar,j A; 



+ 



^Fif ^KL/ \^LI ^FL/ 

^g2A'-4A^^A^^ |^g2A-4A^^A^, _ g2A'-4A^^A„^ 

1 1 \ / A^^ + A,^ 



\^FK ^KlJ K^FK^KL + ^FL^LlJ 

X (e^^'-^^^FK^KL+^FL^Ll) - e'^^'^'^^FK^KL+^FL^Ll)^ ^ (2.22) 

where A^^ is the difference in the invariant-mass of states A and B. 

Next we want to be able to solve for the reduced interaction order by order in 

the coupling, g/<^. However the reduced interaction contains operators and couplings. 
We do not derive this here as it is done in earher work. Simply stated, we use Eqns. and 



(2.13) and match powers of the interaction. 
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This means the change in the reduced interaction when lowering the cutoff is due to 
the renormahzation of the operators and the couphng. If we want to see order by 
order in the couphng how the reduced interaction changes as the cutoff is lowered 
due to the renormahzation of the operators we need to remove the change due to the 
renormahzation of the coupling. This is done next. 

Expand the reduced interaction and 6V in powers of and g^', respectively: 

oo 

v{A) = T^aiv^H^), 
t=i 

oo 

SV = T.9i'SV^'\ (2.23) 



t=2 



where ^*^*^(A) is the 0{gjj contribution to the reduced interaction and 5V^^^ is the 
0{g\,) contribution to 5V . The expansion of 5V starts at second order because the 
first order reduced interaction is the canonical interaction which is unchanged by the 
cutoff. We can now expand Eq. (|2.19| ) in powers of gi^ and giyi: 



t=l t=l t=2 



This equation contains the coupling at different cutoffs, so we expand g\ in powers 
of gK'- 

00 

^7A = ^?A' + E^?1'C.(A,A')- (2.25) 

s=3 

We can use Eq. ( p.25|) to determine g\ raised to the power t > 1 in powers of g\. 

oo 

^7i = ^?l'+E^?i''^M(A,A'), (2.26) 



s=2 



where the Bt^s^s can be determined in terms of the C's by raising Eq. (p. 25 ) to the t 



power. We now substitute Eq. ( p.26| ) into Eq. ( p.24|) and match powers of gA' which 



yields: 

r-l 

V^'\m, A) - V^'^Xm, A') = 6V^'\m, A, A') - ^ Br-s,sV^"~"\m, A), (2.27) 

s=2 

where the -B^-s.s's are functions of m, A, and A' that contain information on the scale 
dependence of the coupling. Since the scale dependence of the reduced interaction 
comes from and the y*^^\m, A)'s [See Eq. (|2.18|) ], Eq. (|2.27| ) simply states that if 



we subtract from 5V^'^\m, A, A') the contribution due to the scale dependence of the 
coupling, then we are left with the contribution due to the scale dependence of the 
V^W(m,A)'s. 

If there is a part of V^'^\m, A) that is independent of the cutoff, it will cancel 
on the left-hand-side of Eq. ( p.27| ). For this reason, we split V^'^\m,A) into a part 



that depends on the cutoff, VQY){m, A), and a part that is independent of the cutoff, 

F("\m, A) = vS>{m, A) + V^\\m). (2.28) 

We must solve for both VQ^{m, A) and V^j^m) to find the IMO. The recursion relations 
for I^Sm, A) and V^\\ m) are given in Sections |2.5.2| and p.5.3| , respectively. 



2.5 Addition of Particle Masses 

This method of renormalization has been generalized to include particle 
masses and was illustrated using massive cjr' theory in 5+1 dimensions . This theory is 
asymptotically free and its diagrammatic structure is similar to QCD, which make it 
a good perturbative development ground. It is straightforward to extend the method 
for massless theories developed in Ref. to calculate QCD quantities for which 



particle masses are unimportant, such as the low-lying glueball spectrum [^, p7 |. 
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In this section, we show how to incorporate particle masses non-perturbatively as a 
necessary step toward a treatment of full QCD. 

In our renormalized scalar theory m is the physical particle mass. In a confining 
theory m is the particle mass in the zero-coupling limit, although other definitions 
are possible. Since the mass is being treated non-perturbatively, it must be included 
in the free part of Ai'^{g^,m, A) in Eq. ( p.l6|) . This alters the unitary transformation 
and leads to fundamental changes in the renormalization procedure. 

The changes in the procedure are discussed in the next three subsections. The 



redefinition of the coupling (Sec. p.5.1| ) is straightforward. In Sections |2.5.2| and 
2.5.3| , we present the expressions for the matrix elements of V^^\m, A) and VQj\m), 



respectively. We also qualitatively discuss the additional steps that are required to 



interpret and use them in a massive theory. In Section |2.5.5| we present the changes 
to the procedure when generalizing from massive 0^ to full QCD. 

2.5.1 Coupling 

The coupling gives the strength of the three-point canonical interaction in the 
theory, thus the canonical definition of the coupling is 
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^203 1 -ML 



l) IP2=P3- 



(2.29) 



647rV5(5) {p,-p,-p,)] 'exp (^9^j (0203|-M'(g^, m, A)|0i) 
GAir^ptd'-^^ (pi -p2- Pa) 



IP2=P3 



In the massive theory, we choose 

[ ' " ^~\0203|l^(^^,m,A)|0i)|p,=p3, (2.30) 

which differs from the definition in the massless theory by the factor exp (9^^. This 
choice of coupling cancels the added mass dependence in the regulator [Eq. (|2.17| )] 
and allows us to closely follow the formalism developed in the massless theory. 
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2.5.2 Cutoff-Dependent Contributions to V^'^Xm, A) 

Momentum conservation implies that any matrix element of V^'^\m, A) contains a 
sum of terms, each with a unique product of momentum-conserving delta functions. 
Assuming that approximate transverse locality is maintained, the coefficient of each 
product of delta functions can be written as an expansion in powers of transverse 
momenta. In massive (p^ theory, we can also make a generalized expansion in powers 
and logarithms of m. For a theory in 5+1 dimensions, the Hamiltonian has dimensions 
[A^], where A is the cutoff on change in invariant mass. For the dimensions to work out 
properly, the scalar fields must have a dimension of [A^] . Thus, the scale dependence 
of any term in this expansion has the form 



where Ni^t is the total number of particles in the final and initial states that partic- 
ipate in the interaction. Also a, (5, and 7 are non-negative integers. For simplicity 
we display one component of transverse momentum, p^; however, the general form 
includes a product of all transverse components from all particles. In principle, the 



the term is independent of the cutoff and is referred to as a "cutoff-independent 
contribution. These contributions are discussed in the next subsection. 




(2.31) 




6 - 2N,^, - a - 7 = 0, 



(2.32) 
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The expression for a matrix element of V^D^m, A) is derived from Eq. ( p.27 ): 



{F\vS>im,A)\I) = [(F|5VM(m,A,A')|/) 



r-l 



(2.33) 



s=2 



- A terms 



"A terms" means the terms in the momentum and mass expansion that contain A' 
are to be removed from the expression in brackets. In terms that depend on positive 
powers of A', we do this by letting A' 0, and in terms that depend on negative 
powers of A', we let A' — > oo. 

2.5.3 CutofF-Independent Contributions to V^'^X'm, A) 

Considering the condition in Eq. (|2.32|) , only two-point and three-point interac- 
tions can have cutoff-independent contributions. The lowest-order cutoff-independent 
three-point interaction is VQj\m) and has not been explicitly computed in the mass- 
less and massive theories. However, Vcj\m) is the lowest-order cutoff-independent 
two-point interaction and must be calculated before anything is calculated to third 
order. 

Due to boost invariance, (m) must be independent of the interacting particles 
transverse momentum, implying 7 can only be zero.Q This means, for example, 
the cutoff-independent part of a self-energy contribution will be proportional to m^. 
Thus, to isolate the cutoff-independent part of a matrix element we must expand it 
in powers of ^ and in powers of log (x)' Then the term that is independent of the 
cutoff obeys the relation: 



"'^''In light-front coordinates, a transverse boost shifts all transverse momenta. This means to ensure 
boost invariance VJjj (m) must be independent of the interacting particle's transverse momentum. 



6 - 2Nint -a = Q. 



(2.34) 
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The matrix elements of VQi\m) are divided into 2-point and 3-point contributions, 
and are given by the expression 

{F\vi\\m)\I) = 

'{F\6V^'+^\m,A,A')\I) 

r,2 

-1 3— point 

s=3 J m^p^ term 

-^[(F|5r(^+V,A,A')|/) 

-Dr,2 '- 

-1 2-point 

-Y.Br+2-s,smV^^''''~'\rnm . (2.35) 



B. 



3 J m"^ term 

Here, ^^mPp^ term" and "m^ term" mean expand the term in brackets in powers of 
external transverse momenta and in powers and logs of m, and keep only the term 
that is proportional to mPp^ or m?, respectively. The removal of A and A' dependence 
is guaranteed by construction. 

Initially Eq. ( p.35| ) looks useless because V^i^m) depends on VQi^^\m) [which 
is inside an integral in suggesting the theory must be solved to all orders 

simultaneously. However, contributions to the reduced interaction from three-point 
interactions can only appear at odd orders, and contributions from two-point inter- 
actions can appear only at even orders. Thus, in the massless theory, this apparent 
problem does not manifest itself because there are no cutoff-independent two-point 
interactions. In the massive theory, although there are cutoff-independent two-point 

(2) (3") 

interactions, it is possible to solve for VQi\m) and V^i (?ti) simultaneously, without 
considering higher orders. This even-order /odd-order pattern can be extended to all 
orders. 
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Including self-energy contributions, the theory we want to describe contains par- 
ticles of mass m. We can simplify the problem by using this fact instead of using 
Eq. ( |2.35| ) to solve for the even-order Vcj\m)'s. We do this by forcing the com- 
pletely disconnected parts of the forward T-matrix elements to be zero. (This part 
of a T-matrix element contains initial and final states that have the same number of 
particles, n, and n momentum-conserving delta functions.) This fixes the even-order 
Vcj\m)'s since they only involve interactions on single particle lines. This allows us 
to calculate (m) independently of V^j (m). This extra condition can be used to 
fix all even-order VQj\m)'s. 

2.5.4 Results 

The coupling in massive 0^ theory runs at third order. We can compare the 
coupling at two different scales, A and A': 

oo 

^/^ = ^/^,+E<,C',(m,A,A'). (2.36) 

s=3 

We can determine how the coupling runs at third order by solving for C3(m, A,A') 
(which is proportional to the matrix element (0203|5^*-'^\'^, A, A')|0i)|p2=p3). Figure 
|2.2| shows how C3(m, A, A') depends on the mass. The running of the coupling is 
exponentially damped as the mass grows since the cutoff inhibits production of inter- 
mediate particles. The difference between the values of the running coupling at two 
different scales increases as the two scales are separated. This is shown by the larger 
magnitude of C3(m, A, A') as the separation between A and A' grows. 

Determining VQj\m) requires a fifth-order calculation and is not attempted. How- 
ever, calculating the matrix element (0203| VqdI"^' ^)I0i) gives the relative sizes of the 
non-canonical interactions and the canonical interaction. Their relative magnitudes 
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Figure 2.2: The third-order coefficient of the running couphng as a function of the 
particle mass. Curves for various upper cutoffs with fixed lower cutoff show the 
coupling is exponentially damped with increasing mass. 
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Figure 2.3: The matrix element of the non-canonical part of the invariant-mass op- 
erator for (pi —>■ 0203 versus the magnitude of the relative transverse momentum of 
particles 2 and 3 in the center-of-momentum frame, y is the longitudinal momentum 
fraction carried by particle 2. 
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are similar to those in Ref. [25|, suggesting that an expansion of the reduced interac- 



tion in powers of the running couphng is vahd through third order. 



Figure shows how the non-canonical part of the matrix element of the invariant- 
mass operator for the interaction 0i 0203 depends on the magnitude of the rela- 
tive transverse momentum in the center-of-momentum frame. Increasing the relative 
transverse momentum in the center-of-momentum frame increases the free mass of 
the system. 

2.5.5 Generalizing to Full QCD 

As in the 0^ theory, momentum conservation implies that any matrix element of 
V^^Xm, A) in full QCD contains terms that can be written as an expansion in powers 
of m and powers of logarithms of m. However, in QCD we work in 3+1 dimensions, 
so the Hamiltonian has dimension [A^] and the particle fields have dimension [A]. So 
for QCD, the scale dependence of any term in this expansion has the form 



/ m 
VA 



m 
log J 



a) ' 



This is the analogous to Eq. |2.31| in 0^ theory. Thus if /3 = and 



(2.37) 



4 - iVi,t - a - 7 = 0, (2.38) 

the term is independent of the cutoff. Thus for QCD there are cutoff independent 
interactions for two, three and four interacting particles. However, since we are 
approximating mesons as a color singlet qq pair, there can only be either two or four 
interacting particles. So the expression for the cutoff independent interactions is given 
by: 
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{F\vS>im)\I) = 

Br,2 



^[(F|5\/(^+V,A,A0|/) 



-1 4-point 

•5=3 J m^pj term 

+ ^\{F\6V^^^'\m,A,A')\I) 

-Dr,2 

r+1 -| 2-point 

-Y.Br+,_,,,{F\V^^^'-%m)\I) . (2.39) 

*=3 J term 

The expression for {F\ VQ^{m, A) |/) is the same in full QCD as it is in (p^ theory. 
2.6 Mass- Squared to Second Order 

We want to find which matrix elements we will need to calculate to find: 

{F\M'{g„m,A)\I) = MUF\I) + e-^{F\V{g,,m,A)\I) (2.40) 



to second order in g\. Thus, from Section 2.4 we must determine 



{F\V^'\g,,m,A)\I) = {F\v!^{m, A)\I) + {F\vg^ im)\I), (2.41) 

where the second-order reduced interaction is divided into the part that depends on 
the cutoff (CD) and the part that is independent of the cutoff (CI)!^. 

There are only three types of diagrams we must consider when calculating the 
reduced interaction to 0{g'^). The self-energy (SE) is a one-body operator that acts 



on both quarks and antiquarks. These diagrams are illustrated in Figure |2.4|. There 



are two types of two-body exchange diagrams. The first is the instantaneous gluon 



exchange (IN) and is shown in Figure |2.5| . The second type of exchange is a single 
"See Eq. (|3|) for the rule to deter mine cutoff dependence or independence. 
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Figure 2.4: The quark and antiquark self-energy diagrams. 



2 



Figure 2.5: The instantaneous exchange diagram. 



gluon exchange (EX) from quark to antiquark or vice- versa. These exchange diagrams 



are shown in Figure 2.6 



The self-energy (SE) interaction has a cutoff-dependent and cutoff-independent 
part, so we split it into these two parts: 



{F\V^{m,A)\I) = {F\V^Z.i^,A)\I) + {F\V^^>Jm)\I). 



(2) 



(2.42) 



The second-order cutoff-dependent reduced interaction can be written: 



{F\V^'>,{m,A)\I) = {F\6V^'^\I) 



A terms • 



(2.43) 





Figure 2.6: The gluon exchange diagrams. 
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For two-particle final and initial states, the matrix elements of {F \ 5V^'^^ \I) are given 
by: 

= {F\V^l^{A)\I)Tt''''\F,KsE,I)\Atern.s 

+ (F|Fg(A)|/)Tf'^')(F,irEX,/)|Aterms, (2.44) 

where Kse, and K^x, are the intermediate states associated with the quark self-energy 
and gluon exchange, respectively. 

The cutoff-independent contribution comes from the instantaneous gluon exchange 
and the cutoff-independent part of the self-energy. The instantaneous exchange is 
cutoff-independent since there are four interacting particles, and there are no terms 
proportional to any nonzero power of the quark mass or external transverse momen- 
tum. Finally we write the cutoff-independent part of the reduced interaction as: 

{F\vg\m)\I) = 

{F\vl^\m)\I) + {F\V^l^im)\I)Tt^'''\F,KsE,I)\Aterms. (2.45) 

2.7 Combining the Interactions 

This section describes how the interactions are combined to explicitly cancel di- 
vergences in the context of the meson calculation (quark- ant iquark states), however 
the general discussion is true for the glueball calculation (glue-glue states) with the 
quark and antiquark replaced with gluons. We list how the interactions are divided 



and recombined which is simply the overview of Perry's work in . 

The self-energy, gluon exchange and instantaneous gluon exchange diagrams pro- 
duce divergences when the longitudinal momentum of the exchanged gluon vanishes. 
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These 'infrared' divergences are an artifact of the hght-front quantization and there- 
fore should cancel. Although the divergences cancel, Perry uses the QCD Hamiltonian 
to 0{a) to show how a logarithmic confining mechanism arises after the cancellation. 

We regulate the divergences with a cutoff on longitudinal momentum. However, 
after the matrix elements are combined, as described below, this cutoff can be taken 
to zero and all matrix elements in the renormalized Hamiltonian are finite. The com- 
bination of matrix elements cancels the apparent divergences so no renormalization 
is needed to remove divergences from small longitudinal momentum. 

The divergent parts of the exchange interaction and the self-energy are cancelled 
by different parts of the instantaneous exchange. We can make this cancellation 
explicit by dividing the instantaneous interaction into two parts, one above the cutoff 
and one below the cutoff: 

{q3q4\M\A)\qiq2)m = {q3q4\M\A)\qiq2)tf, + {q3q4\M\A)\qiq2)l„ (2.46) 

where 

{q3q4\M'iA)\qm)tN = (l - e^^"^--^-) (M4|A1'(A)|M2)in, 
{q3q4\M\A)\q^q2)f^ = e2^-'^--^-(?3?4|A^'(A)ki?2)iN. (2.47) 

Next break up the self-energy and exchange terms into finite and divergent pieces: 

(g3g4|M^(A)|gi^2)sE = (M4|A^^(A)|?ig2)sE+(?3?4|M^(A)|?ig2)|E> 

(M4|A^'(A)|gig2)Ex = (M4|Al'(A)|gig2)gx + (M4|A^'(A)|M2)ex- (2-48) 

We combine the divergent part of the exchange interaction and the instantaneous 
interaction above the cutoff into the finite expression: 

(g3?4|Al^(A)|gi^2)iN+EX = (g3?4|Al^(A)|gig2)ix + {Q3qi\M'^{A)\qiq2)j^. (2.49) 
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We also combine the instantaneous interaction below the cutoff with the divergent 
part of the self-energy because the divergences in the two terms cancel. Thus we 
have, 

(M4|At'(A)kig2)f/ = (M4|M'(A)ki?2>fN + (?3g4|A^'(A)ki?2)FE, (2.50) 
leading to 

{qm\M'^{^)\m2) = (M4|-M^(A)kig2)KE + {q3q4\M^{A-)\qiq2)lE 

+(g3g4|A1^(A)|gig2)Ex + (g3g4|A^^(A)|gig2)iN+Ex 

+ {q3q4\M\A)\qiq2)f/ (2.51) 
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CHAPTER 3 



Basis Functions: B-Splines 



In this Chapter we give some background on the variational method as discussed 
in most quantum mechanical texts. We then discuss the approximation of functions 
using basis functions. Next we introduce B-splines, and motivate their use for a set of 
efficient basis functions. We conclude the introduction to B-splines with some simple 
examples. The first example is simple function approximation and the last two are 
the one and two-dimensional harmonic oscillator. 

This Chapter introduces the basis functions we use to represent the longitudinal- 
and transverse-momentum degrees of freedom in our bound-state calculations. How- 
ever, a good understanding of B-splines is not required to follow the meson calculation. 



If a detailed understanding of the B-splines is not needed, Section |3]^ can be skipped. 
The most important facts about B-splines for this calculation are that they have a 
finite nonzero range and that each B-spline only has a spatial overlap with a limited 
number of other B-splines. 

3.1 Eigenstate Approximation 

In this section we show that when using a particular set of basis functions to 
approximate eigenstates, the lowest eigenvalue of the Hamiltonian in the approximate 
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basis is never less than the real ground state of the system. If more basis functions 
are used, giving a better approximation of the real ground state, the eigenvalue of 
the lowest approximate state should decrease and converge to the real ground state 
eigenvalue. Finally, as the number of basis functions increases the lowest eigenvalue 
should converge to the ground state eigenvalue. The following discussion is derived 
from p4 |. 



The expectation value of the operator H in the state is given by: 

where Eq is the smallest eigenvalue of H. If is expanded in the set of eigenfunctions 
of H, 

\^)=Y.Cn\ct>n), (3.2) 
n 

the expectation value can be written: 

= \cn?E^ >EoJ2 (3.3) 

n n 

with 

m)=j:\^n\'. (3.4) 

n 

The only way for the equality in Eq. (|3.3| ) to be true is if all of the c„'s are zero except 
for Co, indicating |0o) is the ground state. We want to approximate each eigenstate 
of H with a finite set of basis functions \Bi): 



N 



f 



.) = E«ri5.), (3.5) 



Li 

1 = 1 



where a\ is the coefficient of the i basis function when approximating the n 
eigenstate and Nf is the number of basis functions. 
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If one more basis function is added to the set {Bi}, will be a better (at worst 
the same) approximation to the eigenstate. If this extra function could produce a 
worse approximation, its coefficient, a-"^ will be zero. Thus, as more basis functions 
are used to approximate the ground state, the lowest eigenvalue will converge to the 
ground state eigenvalue. However, if the original set of basis functions is altered 
when adding the new function, the lowest eigenvalue may not decrease each time a 
new function is added. 



The success of a Hamiltonian approach will depend on the choice of basis functions. 
The basis functions are used to approximate the real state of the system, so if these 
functions are very different from the real states, it will take a large number of functions 
to approximate the real state. Thus, the convergence of the approximate state to the 
real state is slow if a poor basis is used. A large number of basis functions leads to a 
Hamiltonian with many matrix element s.p^ In this calculation, the matrix elements 
are determined by numerically calculating five-dimensional integrals, which is cpu 
intensive. Therefore, a choice of basis functions that limits the number of integrals 
that need to be calculated is important. 

We can approximate the function f{x) by using a finite set of basis functions: 



""^^The number of matrix elements is proportional to the number of basis functions, for each degree 
of freedom, squared. 



3.2 Motivation for B-Splines 




(3.6) 
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Since f{x) and gi{x) are known, the coefficients are determined by multiplying 
both sides of Eq. ( |3.6| ) by gj{x), integrating over x, and solving: 



dxf{x)gj{x) =aij dxgi{x)gj{x), (3.7) 

where the sum over i is implied. If the gi{x) are orthogonal, then the right hand side 
is nonzero only for i = j. 

To determine the matrix elements of the operator, 0{x), in the approximate basis 
[using the (?j(x)'s], it is necessary to compute integrals that look like: 

dx 0{x)gj{x)gi{x). (3.8) 

These integrals are generally nonzero even if z 7^ j. The number of integrals that 
need to be determined can be reduced by choosing a set of basis functions that are 
non-zero over different ranges of x. 

The basis functions known as B-splines (Basis Splines) have the property that they 
are non-zero only over a finite range of their argument, and this range is different for 
each spline in the set. However, they are non-orthogonal, which means the right hand 
side of Eq. ( p.7|) has non-zero off-diagonal terms. 

3.3 Introduction to B-Splines 

In this section we introduce B-spline functions and outline their derivation by 



Niirnberger and de Boor We begin with a discussion of the "knots" or 

"control points" that determine each spline's shape and then discuss a few of their 
basic properties. Then we state the recurrence relation for the B-splines as well as 
their polynomial generators, discussing subtleties that should be understood. Note 



all of the equations and notation are introduced in reference When using these 
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two references, there are subtle differences in the notation that can cause confusion 
when comparing derivations. 

Since the goal of this section is to give a basic understanding of B-splines, we write 
the definition of a B-spline and then return to the derivation. The i'^^ B-spline can 
be written: 

i+m+l 

Brii)= E aj{t-x,re{t-x,)e{x,^rn+i-t). (3.9) 
j=i 

t is the argument of the B-spline and Xj is the j'*'* knot, aj is a numerical coefficient 
for which we must solve. This B-spline is made of m*'^ order polynomials, thus the 
index m gives the order of the B-spline. The index of the B-spline, i, is associated 



with the knots that are discussed in Section p.3.1| . The range of i is — m < i < k 
where k controls how many B-splines make up the basis, giving m + k + 1 B-spline 
functions in the set. 

3.3.1 Knots 

The knots in a given knot sequence are labeled: 

X^rn < . . . < X^i < a = Xq < Xi < . . . < Xk < Xfc+l = b < ... < Xk+m+1, (3.10) 

where a and b define the range of t for which the B-splines will form a basis, m is 
the order of the B-spline and k allows us to choose how many individual splines we 
use to form our basis. The knots do not need to be equally spaced but must be non- 
decreasing.]^ It is possible to place multiple knots at one point, but we do not discuss 
this until section ^.3.4|, because it complicates the following recurrence relations. 



^•^If there are regions of phase space that are known a priori to be more important (the functions 
being approximated may have more structure in these regions) than other regions, clustering knots 
in the important region produces a set of B-sphnes that will approximate the real functions with 
fewer B-splines. 
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3.3.2 Basic Properties 

One of the benefits of using B-splines for basis states is tliat they are non-negative 
and have a spatial overlap with a limited number of other B-splines (the number of 
overlapping splines depends on the order). The i*'' B- spline of order m, BY^{t) is 
positive in the range [xj, Xj+m+i], and zero outside. Explicitly: 

Bl^{t) =0, t> Xi+m+i or t<Xi (3.11) 
BTit) > 0, Xi<t< Xi+m+i- (3.12) 

For a given knot sequence, the set of m + + 1 B-splines {B^^, . . . , B^} forms a 
basis (they are linearly independent) on [a,b]. 

The coefficients aj in Eq. p.9| ) can be found by solving the linear system of 
equations given by: 

0^x^ = 0, r = 0,...,m (3.13) 

j=i 

i+m+l 

J2 a,x-+i = (-l)-+i(m + l), (3.14) 

j=i 
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3.3.3 Normalized B-Splines and the Recurrence Relation 

The normalized B-splines are defined such that the sum of all B-splines at a given 
point is 1, 

k 

E Nnt) = 1. (3.15) 

i=—m 

The normalized and unnormalized B-splines are related by: 

Nnt) = ^-(x,+^+i - x,)Brit). (3.16) 
m + 1 
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The recurrence relation for the normahzed B-sphnes is: 

iV-(t) = -l—^Nr-\t) + ^!±!!^±LZ^iV^-i(t). (3.17) 

Finally, the n^^ derivative (designated by the "(ra)" superscript) of a normalized B- 
spline is given by the recurrence relation: 

m = (ivr*)'""' w — Ki')'""' (*). (3.18) 

3.3.4 Degenerate Knots and the Recurrence Relation 

It can be useful to place multiple knots at the same point (degenerate knots). 
Although there can be degenerate knots in the region [a, b] we choose to only use 
degenerate knots at a and b. This is because degenerate knots create discontinuous 
derivatives at the degenerate knot. If these discontinuities occur at the boundaries, 
then one can still safely take derivatives in the region of interest, and take the deriva- 
tive at a or 6 to be the limit as the point is approached from the right or left, 
respectively. The recursion relation given in Eq. ( |3.1?| ) must be carefully applied 



with degenerate knots because the denominators can be zero. However, each term on 
the right-hand side of Eq. ( ^.17] ) is finite because in the limit a denominator becomes 
zero, the product of the denominator with the normalized B-spline is finite. 

It should be noted that using splines with degenerate knots can lead to a non- 
Hermitian Hamiltonian. Consider the nonrelativistic kinetic energy term in position 
space: 

H,, = I' dxBi{x)-^Bj{x). (3.19) 
If we integrate by parts, this becomes 

Hij = Bi[x}—-—\^=b - Bj[x}——\^=a - / —r T— . (3.20) 



dx dx Ja dx dx 
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It is obvious the last term is unchanged if we let i ^ j. 

The first {i = —m) and last {i = k) B-splines do not go to zero at a and b, 
respectively, 

B-M ^ 0, Bkib) ^ 0. (3.21) 

This fact prevents H from being hermitian. Consider the first (left-most) spline 



\i = —m): 



.dBAx) ^ . .dB^rr,(x) 



dx dx 
Since j ^ —m we can rewrite Eq. ( |3.22|) as: 



(3.22) 



dBj{a) 

H~mj — Hj-m = -B-.m(a.) — — , (3.23) 



where the expression is evaluated in the limit x —>■ a. The only j for which the 
Hamiltonian is not symmetric under interchange of indices is j = —m+1. For all 
others the derivative at the boundary is zero. So we find for both the first (left-most) 
and last (right-most) splines: 

H-m-m+l — H^m+l-m 7^ H^^k-l — -f^fc-l,A; 7^ 0. (3.24) 

Thus, the only way we can maintain Hermiticity in position space is to discard the 
two splines -B_m and Bk- The problem can also be avoided by working in momentum 
space which avoids second derivatives, the source of the problem. 

3.3.5 B-spIine Polynomial Generators 



Using the recurrence relations in Section |3.3.3| is straight-forward analytically and 



numerically. We can also evaluate the B-splines using the polynomial generator: 

P^i) = E - XjY, t e [x,,a;,+i]; (3.25) 



r=0 ' • 



40 



where, 



.M(t)= Y: «r^r^W, (3-26) 



i=—m+r 

and 

f if r = 

ar = { , .c-D-aC-i) . (3.27) 

I m + l-r^^ ifr>0 ^ ^ 

Although it is easier to understand the recurrence relations, if we want to repeatedly 
evaluate a B-spline it is faster to use the polynomial generator in Eq. ( p.25D . To speed 



up the calculation further, all the a\ s can be calculated in advance. However, these 



equations, which were taken from are used to represent a spline which is a linear 
combination of B-splines: 

A; 

sit) = Y: a^Nnt). (3.28) 

i=—m 

But we are interested in representing only an individual B-spline. Thus if we want to 
represent the B-spline, then we just set all of the for i 7^ / to and aj = 1. So 
we can write the J*'^ B-sphne on the interval [xj,Xj+i]: 



m \ k 



.^/W = E3 E a^iV— (x,) (^-x,)^ te[x,,x,+,]. (3.29) 

r=0 ' I i=— m+r I 

3.4 Simple B-Spline Applications 
3.4.1 Function Approximation 

A simple problem that illustrates the use of B-splines is function approximation. 
We can approximate a function f{x) as: 

f{x)^J2^Mx). (3.30) 
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To solve for the Oj's multiply both sides by Bj{x) and integrate over x: 

rb rb 

I dxBj{x)f{x) = y^a^ / dxBj(x)Bi{x). (3.31) 

J a ^ J a 

Note the range of integration is [a, h] since that is the range over which we are going 
to approximate the function and the range over which the B-splines are defined. This 
is just a matrix problem of the form: 

a, = 07//,, (3.32) 

where 

rb 

Oji = / dxBj{x)Bi{x), (3.33) 

and 

fb 

fj = / dxBj{x)f{x). (3.34) 

J a 



Figure |3.1| shows the eight normalized B-splines of order m = 3 with k = A over 
the range [a, b] with a = and b = 5. At any value of t there are only m + 1 = 4 
non-zero splines. We want to know how well the B-splines are going to approximate 
functions that will come up in QCD. In the longitudinal direction, functions of the 
form (x(l — x))'^ where d is positive, are common in light-front QCD. Figure |3.2| shows 
the function (x(l — x))^ and the approximation built from B-splines with m = 3 and 
k = 5. The knot sequence is also uniform with separation |^ where a = and 



6 = 1. Finally, figure 3.3 shows the difference between the exact function and the 



approximation using B-splines. 
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Figure 3.1: The B-splines B'^{t) for i between —m and k, with m = 3, and k = 4. 
The knots are equally spaced with separation |^ where a = and 6 = 5. 



1.0E-03 



[x(1-x)]'*5 and B-spllne Approximation 



—^approximate 
—»— exact 




Figure 3.2: The function (a;(l — x)Y and the approximation using B-splines. The 
B-splines are order m = 3 with k = 5. 
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Error 




Figure 3.3: The difference between the function — x))^ and the approximation 
using B-sphnes. The sphnes are order m = 3 with k = 5. 



3.4.2 The One-Dimensional Harmonic Oscillator 

To test the B-splines as a useful set of basis functions, we use them to solve the 
one-dimensional harmonic oscillator problem. The harmonic oscillator Hamiltonian 



is: 



H 



(f 1 2 2 

— H — muj X . 

2m dx-^ 2 



(3.35) 



If we write out the eigenvalue equation for a state that we approximate with B- 
splines we get: 

H\<P) = E\<P) ^Y.(^iH\B^) = Ej2a,\B,). (3.36) 

i i 

Looking at a particular matrix element and dropping the explicit sum, we get: 



[Bj\H\Bi)ai = {Bj\Bi)Eai, 
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(3.37) 



where {Bj\H\Bi) is the Hamiltonian matrix {H) and {Bj\Bi) is the overlap matrix 
(O). This problem can be solved in its current form, as a generalized eigensystem 
problem 

Hx = EOx, (3.38) 

or we can rewrite it as a simple eigensystem problem 

Q-^Hx = Ex. (3.39) 

Since a B-spline of order m has only m — 1 continuous derivatives, it is necessary 
to use at least third-order B-splines. Although higher order B-splines may speed the 
convergence, they also add to the total number of states and matrix elements. For 
simplicity we use m = 3 and only adjust k.^ Figure |3.4| shows the eigenvalues for the 



ten lowest states in units of hu. Note the correct value of the energy is En = {n+^)huj, 
so the eigenvalues plotted approach the correct values. The total number of splines 
used is m + A; + 1, or for m = 3, k + 4. The knots are equally spaced over the range 
[a, b] with a = —5 and b = 5. a and b are determined by finding the smallest values 
that, when increased, do not change the eigenvalues. 

The harmonic oscillator can be solved analytically. The eigenfunctions are Her- 
mite polynomials. We can see how well the B-splines approximate the solution by 
considering the overlap between the approximate and exact solutions. The overlap is 
defined as: 

Iadx(f){x)H{x) 



(3.40) 



where (f){x) is the approximate solution built from B-splines and H{x) is a Hermite 

polynomial. The overlap approaches 1 as a better approximation is made. Figure 
^^Using higher order B-sphnes does not produce convergent results in noticeably less time. 
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Harmonic Oscillator Eigenvalues: m=3 




2 5 8 11 14 17 20 



Figure 3.4: Energy eigenvalues of the harmonic oscillator in units of huo. The total 
number of states is m + k + 1. 



3751 shows the overlap of the first three approximate eigenstates determined by B- 
spline approximation with the known solution to the harmonic oscillator problem, 
the Hermite polynomials. 

3.4.3 The Two-Dimensional Harmonic Oscillator 

The two-dimensional harmonic oscillator is a non-trivial extension of the one- 
dimensional harmonic oscillator. The B-spline basis used for the two-dimensional 
problem is a product of the one- dimensional B-spline basis B^y = B^® By where B^ 
and By are one- dimensional B-spline basis. The energy for the two-dimensional case 
is simply the sum from each individual direction: 

E{n^,ny) = ]^{^x + ^y) + (^x^^x + nyUy), (3.41) 
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Figure 3.5: Overlap of the harmonic oscillator eigenstates with the Hermite polyno- 
mials. 



where h = 1. 



Figure p.6| shows the eigenvalues for the two-dimensional harmonic oscillator using 
third order B-splines in both the x and y directions. For convenience we use the 
same, equally spaced, knot distribution is used in both directions. The frequency in 
the y direction is the same as in the x direction which gives rise to the particular 
degeneracies shown. 



The approximation of the lowest seven states are plotted in figure p.7| . The ranges 
shown are limited to —3 < x,y < 3 because the harmonic oscillator states fall off 
exponentially. Note also the phase of the wave function is arbitrary (note the overall 



negative sign in the lowest state). Figure shows the eigenvalues for a two- 
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Figure 3.6: Eigenvalues for the two-dimensional harmonic oscillator, with uo^ = = 
1, and states labeled by {rix^ny). There are k+m+l states in each direction. 
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Figure 3.7: The approximation of the seven of the lowest energy wave functions for 
the two-dimensional harmonic oscillator with Uy = uj^ = 1. The labels are : (0,0), 
(1, 0), (0, 1), (1, 1), (2, 0), (0, 2) and (2, 1). There are k+m+1 states in each direction. 
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dimensional harmonic oscillator without rotational symmetry using third order B- 
splines in both the x and y directions. For simplicity we use the same, equally 
spaced, knot distribution is used in both directions. The frequency in the y direction 
is twice that in the x direction which gives rise to the particular degeneracies shown. 
The approximation of the lowest seven states are plotted in figure |3.9| . The ranges 
shown are limited to —3 < x,y < 3 because the harmonic oscillator states fall off 
exponentially. Note also the phase of the wave function is arbitrary (note the overall 
negative sign in the lowest state). 
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Eigenvalues: m=3, w_x=l, w_y=2 




4 7 10 13 16 19 22 

k 



Figure 3.8: Eigenvalues for the two-dimensional anisotropic harmonic oscillator, with 
ujy = 2ujx = 2, and states labeled by {rix^ny). There are k+m+1 states in each 
direction. 
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Figure 3.9: The approximation of the seven of the lowest energy wave functions for 

the two-dimensional harmonic oscillator with Uy = 2uJx = 2. The labels are : (0,0), 
(1, 0), (2, 0), (0, 1), (3, 0), (1, 1) and (2, 1). There are k+m+l states in each direction. 
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CHAPTER 4 



Matrix Elements of the Light-Front QCD Hamiltonian 



In this chapter we quantize the Hamihonian. We begin with the canonical hght- 



front Hamiltonian as a starting point. Splitting the Hamiltonian density into a 
free and interacting part, we have: 



where 



n 



n 



ggg 



n 



gqgg 



gggg 



— Hfree + '^^int 



"^int — "^qqg + "^ggg + '^qqgg + '^qqqq + '^ggggi 



-9 -2^. 



1 



(^) {a.d-+m)^p. 



id+ 



+2 



d+ 



y_ rabc rade { J\^3 



+2 



9+ 



4 
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(4.1) 
(4.2) 



(4.3) 
(4.4) 

(4.5) 
(4.6) 

(4.7) 



where the gluon and quark fields are defined in Appendix 0. 
4.1 Conventions for the States 

We order quark creation and annihilation operators, b and before antiquark 
creation and annihilation operators, d and S . If there is an initial state that contains 
a quark and an anti-quark, then the ordering will be b^S\0). However, since the final 
state will be the complex conjugate of an initial state, the ordering in a final state 
will be {0\db. Similarly, if a state contains two quarks, then we have: \q1q2) = 
and |gig2)"^ = (gi'?2| = (0|fe2&i- 

The results for the matrix elements of the Hamiltonian presented in Sections ff]B 
and are analogous to [B^] . 



4.2 Quantization of the Hamiltonian 

The full Hamiltonian can be written: 

H = J dx~d^X± {Hiree + Hmt) • (4.8) 

We find: 

\2 /^2 I 



i^free = J ^1 \ + (^^^ ) [^l^l + d\d,] [ , (4.9) 



where Di is defined in Eq. ([A.5|) 

The only parts of the interacting Hamiltonian we need to calculate for an 0{g'^) 
calculation are T^qqg and l^qqqq- We label the transverse momentum of a particle k^^\ 
where a is the particle label and {i) is the direction of the transverse momentum 
component. We define: 

= {s,)k^^ + ikf\ (4.10) 
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and substitute the expressions for the gluon, quark and antiquark fields from |A.4 , 
A.18|, and |A.19| into Eqs. (O) and (O) giving: 



bld^as ^5si,52 + ^S2,S3-^ + ^S2,S3"^| 



-alb2di Ss^,s2 



-imS2Ssi^S2^S2,ss I T+ + T+ 



''I 



Sl,S2"S2,S3 



1 



+5^{ki -k2- k3)x 



dlaldi iSsi,s2 ^ 



h.*S'j, l,*S3 

r ^ r _ ^ 

fs2,S3 7 + ^S2,S3 , + 

1^2 "'l 



—ims25. 



Sl,S2'^S2,53 1 / + 



1 



1^2 '^1 



Sl,S2 



7 S3 J,S3 

"T '^S2,S3 7 + "I" "^52,53 7 + 
"-2 "-1 



k 



1171825 si^S2^S2,S3 \ 7 + 



- A;2 + fc3)x 



b\alb2 ( 5, 



i,*S3 

'^3 



Sl,S2 



6 



l,*S3 
1^2 

S2,S3T4~ 



(5, 



*S3 



S2,S3 7 + 



+zms25.i,s2'^.2,s3 I " A^}) 
/ [ A;^^ 1 

+ d2dia-i ( ^Sl,S2 ~ ^S2,S3"p^ ~ 5s2,^S3"p^| 

+ims25si,S2^s2,s3 I 7~F ~ 7~F 

I K2 /C]_ 



(4.11) 



and 
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mm 



'Nlg^TC' I D1D2D3D4 (ci I T„ I C2) (C3 1 T„ I C4) ^Jktktktkt 

h\hyb2hi + d\d\did: 

6163^264 — d\h2did^ 



\^5\k,-k2^k,-k,)^^4^^ 

+5\k, + k2 + k,-k,)^-^^^ 

{kt + ktf 

-S'{k^-k2-ks + k,) ^^-^^'^^3,S4 



-S^{ki-k2 + k3 + ki) 
+S^{ki + k2-ks + ki) 



61 63^462 — dlb^dids 
bidldlds — 6362^4(^1 



{kt - kty 

+(5^(^i + ^2-^3-^4)%V5%f 

[kj + k^Y 

+6,{k, - k2 - k, - k,)^'''''^'''-' 



bid\b2ds + bld2b4^di 



bUlb^da + bldlb2di 



(kt - kty 



6|62&4<^3 — bld^d^di 



.3 qqg Matrix Elements 



We need to find the matrix elements of all qqg vertices. 
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4.3.1 {qaqb\H\g,) 




This gives us: 

{QaQblHlgc) = 

167r'g^2k+k+{ca\T'^h)S'{ka + h- k,) x 



4.3.2 {gc\H\qaqb) 




(4.13) 



{gc\H\qaqb) 



lQ'K'gpktkt{ci,\T'^''\ca)S\ka + k,- k,) X 

/ ( l.*Sc l.*Sc U* 
/ r I |_ r_ ^ |_ r H 



b 



+imSaSs„sA^,s^ \^ + ^])- (4-14) 
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{Qbgc\H\qa) 



4.3.4 {qa\H\q,g,) 



167T^gj2k+k+{ca\T'^\ck)6%ka -h- k,) X 



^ ( l.*Sc U*^<: U*Sc 

r I ^ r H r '^a 



-zmsb(5,,,,,(5,,,,-, I ^ - } ] • (4-15) 




167i'gj2k+k+{ca\T'^\c,)S\ka -k,- k,) X 



' ( i-'*': 

A J -I- A -I- A - -SL. 

V I '^c '^fc '^a . 



- --- 

""'^ 1 fc+ k+ 



. (4.16) 



4.3.5 (ga^7c|^|gfe) 
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{qagc\H\qb) 



lQT:^gpk+k+{ca\T'^\c,)5\ka -h + k,) x 

is 



U*Sc U*Sc h,*S^ 



k+ 



4.3.6 {qa\H\qbgc) 




{qa\H\qbgc) 



16r'gpk+kt{ct\T'^\ca)S%h -ka + k,) x 

(f h^<= h^<: h'"^" 1 

r I r _ r H_ I 



4.4 gg^^ Matrix Elements 



To calculate a gg state to second order, it is only necessary to calculate the in- 
stantaneous exchange matrix element. 

C ► 



There are two different terms for the instantaneous exchange which give us: 
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{qaqb\H\qcqd)iN = -32^77=^ | L>ii^2i^3A(ci| T„|c2)(c3| T„ [04) V^^/^^^J^ X 

S^(k, -k2-h + h) [Sa,A,5,,2Sd,3 + <^a,3<^6,2<^c,4<^d,l] " (4-19) 

The two different sets of delta functions give identical contributions when we use the 
momentum conserving delta function. Thus, we get: 

{qaqb\H\qcqd) = -6Ag\\cd\T f\c,) {Ca\T f\c,) ^ K k+ kj x 

s^k^-k,-k,+h) ^^i''^T:' - im 
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CHAPTER 5 



Matrix Elements of the Renormalized Hamiltonian 



We now determine the matrix elements of our renormalized Hamiltonian in a 



plane-wave basis. Section shows which matrix elements are needed to calculate 



the mass-squared operator to second order. These matrix elements are calculated in 



Sections 5.1 through 5.3 



We find matrix elements of the IMO which is very closely related to the Hamil- 
tonian matrix elements we have already calculated. Specifically, in the center-of- 
momentum frameQ 

= V+n, (5.1) 
where is the total light-front longitudinal momentum. 
5.1 Cutoff Dependent, Non-Canonical Contributions 

The non-canonical contributions to the mass matrix are contained in 
(^3^4|VcD(A)|gig2), or more explicitly in {F\ 5V^*^^) |/) [Eq. (|2.44| )]. For a second order 
calculation, we can write: 

|/) = W{m4\V^'^\K){K\V^'^\q,q,)Tt^'''\F,K,I), (5.2) 

^ K 

""^^Light- front formalism is discussed in Section 2.1 and Appendix 
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where the sum is over a complete set of (intermediate) states and V^^'' is the first- 
order canonical Hamiltonian. The only term in V^^^-* we need contains a qqg vertex 
so \K) can only be a qqg state, and the interaction is either a gluon exchange or a 
quark self-energy. 

In this section, we make the following momentum definitions for various momenta 
in the self-energy and exchange diagrams: 

ki = {xV+,xV + q-l), k2 = {[l-x]V^,[l-x]'P -ql), 
h = {yV+,yV + pl), k^ = {[l-y]V+,[l-y]V~pl), 
/C5 = {zkf , zki + fl) , k-r = {[1 - z]kf , [1 - z]ki ~ fl) , (5.3) 

where fci and k^ are the momentum of the incoming and outgoing quark and anti- 
quark, respectively, fcs is the momentum of the internal quark. Using these definitions, 
we find for the free mass squared of the states: 

MJ|F) = Mj|g3g4) = |^, (5.4) 
MJ|J) = Mf\q,q,) = ^-^ (5.5) 
M}\K) = M^f\q,q,g-j) 



m?{l — z){l — X + xz) + z{l — z)q]_ + (1 — x)f] 



(5.6) 



a;(l — x)Z[l — z) 
5.1.1 Self-Energy 

The quark self-energy contribution to {F\ dV*^"^^ \I) is: 

i j D,DeDr{q,q^\V^p_,^\q,q,gj){q,qe97\VjX^^^ (5.7) 
where the subscript on V^^^ indicates which of the vertices calculated in Section E]3 



are 



used. For the quark self-energy we will need the matrix elements given in Eqs. ( [4.16| ) 
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and ( 4.17| ). However, these are not the complete matrix elements, since there is a 
spectator quark. We use: 



{q5q&97\V^^^\q2qi) = 



f ( U*S7 U*S7 l,*S7 
r I |_ r _ M I r ^5 



/vy A/-|^ A/^ 



1 



and 



(g4g3|^^^^|g5g65'7) = 



' ( h'^'^ h"'' h'^'' 

A J —-L- _|_ A I A - 

\ y '^i '^5 



1 1 



-ims56s3,s5Ss,,s7 { — - — }]■ (5-9) 



Thus the contribution from the quark self-energy diagram can be written: 

{lQn'V+)^62,4 J D.Djk^^/ktktT^hTyih -h- kr)6^{k, - k, + kj) 



7 f 1 r 5 



X ^s„. 



■'S5,S1 ^ 1.+ ^ "^1.^7 7 + ^ "Sl,S7 7+ r 

I '^5 J 

xTf'^')(F,ir,/). (5.10) 
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After integrating over everything but the momentum of the internal quark, and com- 
pleting all sums, this becomes: 

^ / n''''^ . " \ ^ -^^(F. /). (5.11) 

Recall that z = k'^/kf. For the self-energy, the final and initial states have the same 
energy giving: 

A,, = -A,, = _ r-' + 'f(l-^)\ ,5.12) 

XZ[1 — Z) 



and 



Tf K, I) = --^ (e-2A'-^^L _ e-'^-'^Kr) . (5.13) 
Aki ^ ' 



Using the definitions: 



u = 7z r , P= , 

xz{l — Z) xz 

and changing variables from r to u and integrating over u gives: 



-A 



,2 ' ^1,3^2,4 
247r2 



/^|i±^v^[l-erf(7'(x,z))] 

-2v^y (x, z)Ei (1, 7'2(x, z)) }] - [A' ^ A] (5.14) 

for the contribution to {F \ SV^"^^ \I) from the quark self-energy, where 'erf is the error 
function and 'Ei' is the exponential integral. The anti-quark self-energy contribution 
to {F\ (5^(2) |/) is: 

1 P 

- j D5DQDy{q3q4vff_,g\q5qQgr){q5qQg7\Vj%g\qiq2)T2{^^ Kse, I), (5.15) 

where the subscript on V^^^ indicates which of the vertices calculated in Section 
is used. 
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The contribution from the anti-quark self-energy is the same as that from the 
quark self-energy with x ^ 1 — x. So we find the total contribution to 6V from the 
self-energy diagrams is: 

z dz [l + z"^ 



'^l,3'^2,4 . /2 

247r2 



l-z 

+ 



27r [erf (7' (x, z)) + erf (7' (1 - x, z)) - 2] 



7'(a;,;2)Ei(l,7'2(x,z) 

+7'(l-a;,z)Ei (l, ^^(l - x, z))] }] 



- [A' - A] 



(5.16) 



To find (q'3^4| VgE'*(A)|gig2) we have to keep only the A terms and remove the cutoff 
independent part. We remove the cutoff independent part by expanding our answer 
in a generalized series of powers and logs of ^ and remove the part proportional to 



2 

m . 



Keeping only the terms that contain A gives us: 

z dz ( 1 + z'^ 



" 247r2 



l-z 



2tt [erf (7(x, z)) + erf (7(1 - x, z)) - 2] 



[7(x, z)Ei (1, 72(x, z)) + 7(1 - X, z)Ei (l, 7^(1 - x, z))] }] . (5.17) 

Finally, we have to remove the terms that are independent of the cutoff, since they can 
not be determined by calculating 6V. The only part of this expression that contains 
terms proportional to are the error functions. We find the cutoff-independent part 
of this expression to be: 



247r2 



zx{l — x) 



{5.11 



This term is divergent, however the entire term will be cancelled when calculating 
the cutoff independent part of the reduced interaction (Sec. p73| ). Thus the complete 
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expression for {qaQilV^'^ {A)\qiq2) is: 



(2) A M ^ ;^ \ _ A 2 ^1,3'^2,4 

247r2 



(g3g4|K,E(A)|gig2) = -A 



X 



^ . 27r[erf(7(x,z))+erf(7(l-x,z))-2] 
+V8(7(x,2)Ei [l,7^(x, z)" 

+7(l-x,z)Ei [1,72(1 -x,z)])} 



—4m / (iz 



l + z=^ 



zx(l — x) 



(5.19) 



5.1.2 Exchange 



The first exchange diagram contributes to {F \ SV^"^^ 



(5.20) 



where the subscript on V^^^ indicates which of the vertices calculated in Section 
is used. For this exchange diagram, we need to use the matrix elements derived from 
Eqs. (CT ) and (CT) : 



{q^%g7\V^^^qiq2) = I^'k^'V^ J2ktkt {c^\T'''\ci)52fi6\h - + k-j) 



5s. 



''7 



S5,Si^Sl,S7 



1 



(M4|^^'^|gig2^77) = lQ-K^V+j2ktkt{c^\T:'^'\c^)5^^^5\k,-k^ + k^) 



A J ML _ A _ ML _ A Ml I 

* 1 ■''^^ kt itg" J 

+ims45se,s4^s4,s7 



. (5.22) 



The total contribution from this exchange diagram is: 
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(167r3p+)2 J Dj^ktktktktTl{Tl^5\k2 -k^ + h)5\ks -k^ + kj) 
X 5, 



■'S3,Sl 



t.*S7 h,*S7 l,*S7 



kt 



1 1 



kn 



At^ 

f 1 



xTf'^')(F,ir,/). 



(5.23) 



For the first exchange diagram we find: 



a(i) 

'-^KI 



A 



(1) 

FK 



rn?{x — + {xp — yq)"^ 
xy{x - y) 
rn?{x — yY + [(1 — x)p— (1 — y)qf' 
{I - x){l - y){x - y) 



Finally we use 



C7 U Z 



and find for this exchange graph: 

^ x{l - x)y{l - y) 



x-y 



S\ki + k2-k3- ki) 



X I 5ci,C2'^C3,C4 2'^C2,C4'^Ci,C3 ) T2 



(A,A'){1) 



(5.24) 
(5.25) 



(5.26) 



,S3^S2,S4 ^Si,S4^B 

+im (S25si,s3^s2,s4^+- + ■Si'5si,s3^s2,s4^-+)] , (5-27) 



where the extra parameter on T2 simply specifies it is for the first exchange diagram 
and is defined in terms of the A's in Eq. |5.24| . The 5's are given by: 
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-5 



Sl,S2 



x{l -x)y{l -y) 
1 



+ 



x(l-x) y{l-y) 



x-y 



2y-l , 2a;- 1 2 



-p 

y{l-y) x{l-x) 
p-q{si) p-q{s2) 



xy 



{l-x){l-yy 



B_ 



x-y 



_ gS2 gS2 pS2 

~l~ "Sl .So ~l~ "Sl .So 



(1 — — y) \ x — y ^' ^ X 

xy \ X — y 



p 



1 — X ^ ^ y / 



(5.28) 
(5.29) 
(5.30) 
(5.31) 



x(l - x)y(l - y)' 

For the second exchange diagram, we can use the result of the first diagram and 
switch particles 1 and 2, as well 3 and 4 (See figure This means we will have: 



X 1 — X, 



p -p, 



q -q, 



Sl ^ S2, 



S3 ^ S4, 



(5.32) 



but the color indices do not change. So for the second exchange diagram, we get: 

m^(x — y)^ + {xp — yqf 



A 



(2) 



FK 



,(1) 



A 



(2) 
KI 



xy{x - y) 
m?{x — yY + [(1 — x)p — (1 — y)qf' 



{I - x){l - y){x - y) 
making the contribution from the second diagram: 
^ x{l - x)y{l - y) 



A 



(1) 



(5.33) 
(5.34) 



x-y 



1 



X ( 5ci,C2^C3,C4 2^C2,C45ci,C3 ) T2 



^si,S3^S2,S4B-\ — \- ~\~ TTl 5sj^^s^Ss2,S4^si,S4B 

SiSsi,s3^s2,s4B — h)] • (5.35) 
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The extra parameter on T2 simply specifies it is for tlie second excliange diagram and 

is defined in terms of tlie A's in Eq. ( p.33| ). and B are tlie same 

ones defined for tlie first excliange diagram. Eq. ( p.33| ) shows that 

r^(AA'){i) ^p^ i^, /) = -T^^'^'^^^\f, K, I) = K, I). (5.36) 

There is an implicit cutoff on the exchanged gluon momenta, which we have omitted 
for clarity, that prevents an infrared divergence. Explicitly, in the contribution from 
the first diagram, the cutoff is on \x — y\, and in the second contribution, the cutoff 
is on \y — x\. 

The form of the exchange contribution to the approximate states will look like: 
j d^qd^p j dxdye{\x-y\- e)E{x,q,y,p)Xs^ {y,p), (5.37) 

where E{x, q, y,p) is a function only of momentum. If we split this into two terms, one 
where x > y, the other where x < y, and note that the exchange function is symmetric 
under particle exchange, we can rewrite the exchange integration variables in the term 
where x < y to get: 

[XsiS2 {x, q)Xs3Si {y, P) + Xs2Si (1 - a;, -qjXs^sa (1 - Z/, -P)] ■ (5-38) 

If both of the spin-momentum wave functions are either symmetric or anti-symmetric 
under exchange of particles, then the two terms are the same. If one is symmetric 
and the other anti-symmetric, then the integral is zero. Thus we get: 

2 J d^qd^p j dxdye{x - y)e{x - y - e)E{x,q,y,p}Xs^s2ix,q}Xs,sAy,P)- (5-39) 
Therefore the complete contribution is: 
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-167r^7'+^(x - y)e{x -y- e)5^{ki + k2 - ks - k^) 



1 



x{l - x)y{l - y) 
x-y 



X 



+im (S2(^Sl,S3<^S2,S4-S+- + Sl(^Sl,S3<^S2,S4-S-+)] 



X 



^(1) ^(1) 



2A/-4a(i) A*^' 2A-4A^^'' a 



-4a(1) a(1) 



. (5.40) 



To get the "A terms" we can let A' — > oo. This gives us: 

(g3g4|v^iS(A)kig2) = 

-167r3p+^(a; - y)e{x -y- e)S^{ki ^k^-k^- k^) 

^^C2,cAl,C3 + 5c^,C2^C3,c}j x{l - x)y{l - y) 



X - 



X 



<^Sl,S3^S2,S4-^++ + ^^^S1,S3^S2,S4-^- 



Sl,S3^S2,Si-B-\ + 



X 



,(1) 



.(1) 



(5.41) 



\L^FK ^KI/ 

where B±± — {x — y)B±±. 

Since there are four interacting particles in this diagram, a cutoff-independent 
contribution must be independent of both ^ and It is also obvious that the 
only cutoff dependence in the expression is from the exponential. The 1 subtracted 
from the exponential removes all the cutoff-independent terms, so this is the complete 
expression for {q3q4\\4^}{A)\qiq2) . 
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5.2 V^^^: Instantaneous Diagram 



This is a matrix element of the canonical Hamiltonian, and can easily be derived 
from Eq. ( pq ): 



^ V 3 4 "-1 "'2 



(^3+ - 

-5 



y~\Jx{l - x)y{l - y)' 



(5.42) 



So the complete contribution from the instantaneous diagrams is: 

(?3g4|ViS^(A)|gig2) = {qzqi\Vc\qiq2) 

-647r^P+(5^(A;i + k2 - hi - ki)Jx{l-x)y{l-y) 



y.6{\x — y\ 



-•-3,1 -'-2,4"si,S3W, 



(5.43) 



{x-yf 

We can use the symmetry of the external states to rewrite this as (like in the exchange 
diagram) : 

{qm\vi\miq2) = 



-l2M^V^6\ki + k2-k3- ki)Jx{l - x)y{l - y) 



x6{x — y)6{x — y 



(2) 

5.3 V^Nc- The T-matrix method 



^3,1 ^2A^si,S3^Siy 



S2 



{x - yf 



(5.44) 



In this section we derive the second-order, non-canonical, cutoff-independent part 

(2) 

of the reduced interaction. We show that V^c is the same as the cutoff-independent 
part of the self-energy [Eq. (|5.18|) ] with the opposite sign. This means the sum 
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of the self-energy contribution and the cutoff-independent, non-canonical contribu- 
tions yield the self-energy contribution before the cutoff-independent piece is removed 
[Eq. (1^3)]. 



Two particles that do not interact give a zero forward T-matrix element. We will 
set the "forward part" of the T-matrix to zero to fix the second-order non-canonical 
matrix elements. This means we want a particle of mass m in a state to propagate 
as a particle of mass m. In a light-front theory, we can write the T-matrix to second 
order as: 

T (pi +P2) = ^^int(A) + Hi^tW—-^ r^^int(A), (5.45) 

^ ^ Pi +P2 -h 



where: 



= (5.46) 



and h is the total free energy of the intermediate state. So explicitly, we want: 

= {q3q4\V^'\A)\qiq2) + {qsq^lV^'^ 7 ^^^'^Igi^s). (5.47) 

'P+[pi+P2-h) 

If we write: 

^^^^A) = \4g(A) + V^'^ = \4Sg(A) + + V^'\ (5.48) 



we see that the first term on the RHS is what we calculated in Sec. |5.1| , whereas the 
second term is what we are trying to fix, and the last term contains the instantaneous 
diagrams. However, since we are fixing the forward T-matrix, we only want to in- 
clude diagrams that have the same delta-function structure as a self-energy. So more 
explicitly, what we want is: 

= {qsq^lVsE {^)\qiq2) + {q3q4\vic\qiq2) 

+(^3541^^'^— 7 ^ -y^^'^|gig2)|sE, (5.49) 

V+ [Pi +P2 -hj 
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where (q'3^4|Vnc |9i^2) is the second-order, non-canonical, cutoff-independent reduced 
interaction. The third term is similar to the first (before the cutoff-independent part 
is removed), except that relating the first to the third will show that: 



FK _iiff 
e A" e 



FK 

e a4 e ^ 



P5 -Pe -P7 



e A* e A" 



e 



A 



KI 



A 



KI 



2 -^2 1 A terms 



(5.50) 



This means the third term will cancel the first term, except for the cutoff independent 
part in Eq. (|5.18|) . This leaves us with: 



= (93^41 V'nc 1^1 + 4^i^^^m^ f dz- ^ 



247r2 



zxil — x) 



(5.51) 



givmg: 



^1,3^2,4^2 \ r Zdz [l+z"^ 



r{2) 



K2)| 



247r2 



2tt [erf (7(x, z)) + erf (7(1 - x, z)) - 2] 



+v^[7(x,z)Ei(l,72(x,z); 

+7(l-x,2;)Ei (1, 7^(1 -x,z))]}]. (5.52) 
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CHAPTER 6 



The Basis for the Expansion of Real Meson States 



In this chapter we define the basis we use to approximate the real meson states 
in Chapter 0. Before we use these basis states to determine the meson spectrum we 



test them by calculating the glueball spectrum in Section |6.8| . The basis we use to 
approximate the glueball states is so closely related to the basis for the meson states 
that we do not discuss it in detail. 

We start by assuming the meson will contain one quark and one anti-quark. In 
this case, an expansion of the meson state in the free basis looks like: 



Al^2(gig2|^^"(P))|gig2), 



(6.1) 



P is the three-momentum of the state, j is the projection of the total spin along the 
3-axis, Di is defined in Eq. ( [A.5|) , and n labels the mass eigenvalue of the state. The 
derivation of the approximate state culminates with the final expression in Eq. ( |6.46| ). 

Since the color part of the wave function can be separated from the spin and 
momentum parts, we write: 



(6.2) 
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where |T„) is the color part of the state and |r-^"(P)) is the momentum/spin part. 
Now we can write: 

(gig2|*^"(P)) = (ci,C2|T„)(pi,si,p2,S2|P"(P)). (6.3) 
6.1 The Color Basis 

We define the color part of the scalar product as: 

M;^,,^ = (ci,C2|T„). (6.4) 
Using this, the most general color state can be written: 

|T„)= E<,c.kc,gc.). (6.5) 

Cl,C2 

If we apply a color rotation to a color singlet state (which will not effect the scalar 
part), we see: 



IT') = y M" y q' ) 

I n/ / y ci,C2 l^ci^ £2/ 



Cl,C2 

|T„). (6.6) 



From 1^ a color-rotated quark looks like: 

= u\qc^) = IZ|gc3)(gc3|f^|gci) = Y.Uc^cMci)- (6.7) 



C3 C3 



U = exp[— i^bF^] is the quark color rotation operator where F'^ = = are 

the SU(3) generators for the quark representation, and the di, are the color rotation 
parameters. The rotation operator for the anti-quark representation is given by U* . 



So now we can rewrite Eq. (|6.6| ) as: 

in)= E Ml^c,Uc,c,Ul.Mc.q,,). (6. 
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E 

C1,C2,C3,C4 



Now, U* = W^, so: 

in>= E Mlc.Uc,cM,jQc,qc,). (6.9) 

Cl,C2,C3,C4 

But — U~^, so if we use M^^^^^ = 5ci,c2 then: 

C1,C3,C4 

= El?c4?cJ. (6.10) 



C4 



Now we normalize the color wavef unction: 



(T;|T,) = E (^c^^^ke^ge^) = iV' E ^4,444,4 = iV'E^c4,c4 = N,N\ (6.11) 

C4,C4 C4,C4 C4 



Therefore, the normalized color wave function is: 



1'^") = E ^Cl,C2kcife)- (6.12) 



JN 

6.2 Momentum Conservation and Plane Wave Normaliza- 
tion 

Since the ^'s are mass eigenstates, and since the mass operator conserves momen- 
tum, {qiq2\^^^{P)) must be proportional to a momentum-conserving delta function. 
We write: 



1 1 



lQnH^^\P - p,- P2)^==^eA,,c^Q^^''^^ptpt^ilsM^P2). (6-13) 

V 167r'* yl^c 

where ^i^^s2^Pi-iP2) is the spin/momentum wavefunction and 

e, = e{x-e)e(i-x-e), (6.14) 



76 



which is used to cut off infrared divergences. Now we can rewrite Eq. ( p-lj) as: 
11^/" <Pq dx 



E 



where we have used the Jacobi variables: 



-eA.,c.^ll,,{P,x,q)\q,q,), (6.15) 



Pi 



(xP^, xP± + q 



P2 



[[l-x]P+,[l-x]P± + q 



(6.16) 



and integrated over p2. 6^ is used to regulate the light-front infrared divergences. 
These divergences will be explicitly cancelled allowing us to let e ^ before the 
Hamiltonian matrix is diagonalized. 

Next, we want the state to have a boost-invariant plane- wave normalization: 



^^■'"'(P')l^^"(^)) = 167r3p+53^P - P')<^w'^nx 

We can explicitly write out the left-hand side as: 

11 ^ f (fqdxO, d^q" dx'e'^ 



(6.17) 



167r3 Nr 



E 



^Jx{l — x) \Jx'{l — X') 

$:f^5(P', x\ q> )<^ZJP^ ^' <l){<lW2\<li<l2) 



(6.18) 



where we use: 



{qWi\QiQ2) = -{qW2\QiQ2)- 



(6.19) 



We can write: 



{qW2\Qi'12) 



'pt [P^-Pt) iP-P') 

' P+x(l - x)6^ (P - P') 6{x- x') 6^ {q' - q) 



X <^si ,s'2 4i ,c[ ,4 • 



(6.20) 
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So now if we integrate over q' and x, and complete all the sums except Si and S2, we 
get: 

lQn'P+5\P-P') Y: I d'qdxe,^:{X{P,x,q)^Zs.iP^^^^)- (6-21) 



Sl,S2 



This implies: 



S^,<n' = E / d'qdxe,<l>:{X{P,x,q)<l>ZJP,x,q). (6.22) 

S1,S2 



Since |\&-^'^(P)) is uniquely determined by the equation: 



p-(A) ^^"(P) 



which can be written: 



(6.23) 



M\A) ^^"(P)) = ^•'"(P)) , (6.24) 

'^'i"s2(-^' ^' ^) uniquely determined by A1^(A). Since the free-state matrix elements 
of A^^(A) are independent of the total momentum, so is ^^"^^(P, x, g*). Then the 
normalization condition is: 



ciVdx^,$:^;,(x,g)$i;^,,,(x,g), 



(6.25) 



and the state is given by: 



^i"(P) 



1 1 



,C2 



d'^q dx 
xil — X 



A5c.,cMls,{^.q)\qiq2)- (6.26) 



6.3 Fermions, Charge Conjugation, and Exchange Symmetry 



Charge conjugation is a good quantum number because it is a symmetry of the 
strong interaction. We show |^ that applying the charge conjugation operator, C, 
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gives the eigenvalue ±1, depending on whether or not the spin/momentum wavefunc- 
tion is symmetric or anti-symmetric under particle exchange. 
Using some shorthand, we can write the meson state as: 

l*> = ^hjAuc E / dkdk' ^s,sA k') &I4|0), (6.27) 

S1S2 

where ^sis2^i k') is the spin/momentum part of the wavefunction. /j and q are the 
flavor and color indices, respectively. Next, we act on the state with C, giving: 

= -^/i,/24i,c2 I dkdk' ^sis^ik, k') 44|0) 

S1S2 

= -Sf,jAuc.J2 [dk>dk^s,s,{k',k)b\dl\0), (6.28) 

S1S2 

where the particle indices are switched in the last step. This means that if $ is 
symmetric under particle exchange {si ■(-^ S2, k ■(-^ k') then the state is odd under 
charge conjugation, and if $ is antisymmetric under particle exchange, then the state 
is even under charge conjugation. 

Thus, if we consider symmetric charge conjugation states, then the spin-momentum 
wavefunction must be anti-symmetric under exchange, and an anti-symmetric charge 
conjugation state must have a symmetric spin-momentum wavefunction under ex- 
change. 

6.4 Momentum and Spin Wavefunction Bases 

The Spin-Momentum Wavefunction is given by: 

k^) = Yxr^ix, k±). (6.29) 

9=1 
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The x's are given by: 



xr^ = (6.30) 

where s = —s, and they obey the relation 

T.xrxr = K<i'- (6-31) 

S1S2 

The momentum wave function is expanded in a complete basis: 

fc; — 1 kt oo 

n^ix,k^)= Y: E R^aBi{x)B,ik)AM, (6.32) 

l=—mi+l t=—mt+3 a=—oo 

where m^, kt, mi, and ki are the normal B-spline parameters for the transverse and 
longitudinal basis functions.|3 The Aa{4>) are the basis functions for the transverse- 
angular degree of freedom, the Bi{x) are the basis functions for the longitudinal- 
momentum degree of freedom (Sec. 5.4.1|) , and the Bt{k) are the basis functions for 



transverse- momentum degree of freedom (Sec. |6.4.2|) . 
The transverse-angular basis functions are given by: 

1 



A„(0) = -^e"'^ (6.33) 



and have the normalization: 



^\(f)A:'{<f))Aa{(j))=Sa,a'. (6.34) 



^^The normal range of the B-spline index is from —to to k. However, two longitudinal B-splines 
{I = —mi,ki) and three transverse B-splines (t — —nit.,— mi ~\- 1,—mt + 2) are discarded because 
they produce divergent kinetic energies. 
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6.4.1 Longitudinal Basis Functions 

The longitudinal states are functions of x, the longitudinal momentum fraction 
carried by one particle. The symmetry of the problem allows us to choose functions 
that are symmetric or anti-symmetric under particle exchange {x ^ 1—x). B-splines 
are not symmetric functions about x = | a priori. The longitudinal functions must 
be symmetric and anti-symmetric combinations of the B-splines. We choose knots 
that are symmetric about x = ^, so that pairs of B-splines are also symmetric. For 
example, 

B^-J'\x) = Bi-''\l-x), 

= Btt\l-x). (6.35) 

If there are an even number of B-splines they can all be paired in this manner. If 
there are an odd number, the "middle" spline with index 1 = 1= "^+^+^ _ i is, by 
itself, symmetric about x = |. We choose to always use an even number of splines so 
there are an equal number of symmetric and anti-symmetric functions. Otherwise, we 
would expect the basis to better approximate symmetric states than anti-symmetric 
ones. We will write the symmetric B-splines as: 

Bt:f2^{x) = Bt''\x) + Bt.tUl - ^ < I (6.36) 

and the anti-symmetric ones as: 

Bt.iLjx) = Bt''\x) - Btila -x); l> l (6.37) 

We drop the '(a)symm' notation, and all B-splines are implicitly (anti)symmeterized 
as determined by their index, m and k are the usual indices associated with the splines 
used in Chapter ^ and / is used in place of i to signify it is the index for a longitudinal 
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function. We use degenerate knots such that . . . ,tQ = a and tfc+i, . . . , tfc+m+i = b, 
and we will not use the splines -B_m and Bk because they do not have the correct 



6.4.2 Transverse Basis Functions 

We use t as the index for the transverse basis functions. However, the knots that 
determine each B-spline are designated by tj. Therefore, if a t has a subscript, it 
refers to a knot, otherwise it is a transverse function index. 

The magnitude of the relative transverse momentum in a state can be in the range 
[0, oo). However, to compute integrals numerically we need to restrict this to a finite 
range. We change variables: 



The range of z is restricted to [—1, 1]. This finite range should lend itself well to 
using B-splines, as should the fact that the integrand under this change of variables 
tends to be reasonably smooth. We want a set of basis states that after this change 
of variables produces B-splines. So let us start at the end, and work backwards to 
develop this set of basis states. 

We use degenerate knots such that 



behavior]^ as x — > 0, 1. 



Finally the longitudinal states obey the normalization: 






(6.40) 



17i 



The basis functions must vanish at a; = 0, 1 to ensure finite kinetic energy. 
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and 

^fc+l, • • • , tk+m+l = b. (6.41) 

We are interested in the range z G [— 1, 1]. Choose the knots t_m, ■ ■ ■ ,to = — 1, ti, . . . , tk+i 
1, . . . , tk+m+i, with the corresponding set of B-sphnes {B_rn, ■ ■ ■ , B^}. Now spread 
these states out over the range [0, oo) using the change of variables in Eq. (|6.39| ). 
The knot distribution is changed so that 

t; = ^ - 1. (6.42) 

This smeared set of sphnes will be labeled Bt, t G [—m, k]. The kinetic energy part 
of the Hamiltonian after the change of variables is: 

showing that as limj^^_i, B^'^{y)B^'^{y) must die faster than (1 + yY. This requires 
the use of splines of order m > 3. However we must also consider the power series 
representation of the B-spline, (which is non-zero between to and ti). We write 
it as: 

m 

BT{y)= E + (6.44) 

p=m+t 

where t can be negative, and this form can best be understood by considering Eq. 
( p.25| ). This means that to ensure there is no divergence, the lowest B-spline we can 



keep is B^^_^{y). The transverse basis functions obey the normalization: 

dkk\Btik)\^ = I. (6.45) 





We write the approximate meson state, built from our basis functions as: 

^^"(i')> = ^E^'.rtk,M,j), (6.46) 

qlt 
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where 

xA,_,,_,,(0)|gig2). (6.47) 



Finally, N is used to ensure plane wave normalization (Eq. |6.17|) . 



6.5 Restricting the Spin-Momentum Wavefunction 

If we consider states with different charge conjugation separately, we can re- 
strict the spin-momentum wavefunction to be either symmetric or anti-symmetric 
under exchange. For a negative charge conjugation state, it is necessary for the 
spin/momentum wave function to be symmetric under particle exchange. As de- 
scribed in section |6.4.1| , the longitudinal states, Bi{x), are symmetric under x 1—x 



if / < / and anti-symmetric if I > I (recall I refers to the "middle" spline). If a is even, 
then the angular state is symmetric under exchange, and odd if a is odd.^^ Finally, if 
g = 1, 2, 3 in the spin state it is symmetric under exchange, and it is anti-symmetric 
if g = 4. Thus, for the entire spin/momentum wavefunction to be symmetric we must 
have one of the following conditions: 
For j even: 



q = l,2,4^/>/, 

q = 3-^l<l, (6.48) 

For j odd: 

q = 1,2,4^/ <r, 

""^^By requiring the approximate meson wavefunction to be an eigenstate of rotations about the 
3-axis, AUen |^ proved that a — j — si ~ S2, where j is the projection of the total spin along the 
3-axis. 

84 



q = 3 ^ Z > Z. (6.49) 

Next for positive charge conjugation states, we want the total spin-momentum 
wavefunction to be odd under exchange, giving: 
For j even: 

q = l,2,4^i<r, 

q = (6.50) 

For j odd: 

q = 1,2,4^ Z>r, 

q = s^l<l. (6.51) 

6.6 Meson Overlap Matrix 

Since the B-sphnes are non-orthogonal, we need to calculate the overlap matrix. 
We use this section to carefully derive the form of the overlap matrix and find relations 
that will be useful in the two-dimensional Hamiltonian integrals. 

A free state with spin-state index q, and spin j will be written: 



11 r 



S1S2 



C Sl,S2,Ci,C2 ' 



x^M=B,{k)A^_,,_,,{<t>)\q,q^). (6.52) 

^Jx{l-x) 



Then: 

1 1 



{q',l',t',j'\q,l,t,j) = -— — Yl f d'k^d'k'^dxdx'eA' 

iOTT iVc s-i_,S2,Si, ,S2/ ,ci,C2,Cii ,C2/ 

x5c„c.4„e„X^^^^X;r^'^*_,,_,,(0)^,v_,^,_,^,(0') 

x^£M£d^^B,{k)Mk'){qi'q2'\qiq2) (6.53) 
\/x{l — x)x'{l — x') 
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where 

(gi'?2'kl?2) = (^l,l'(^2,2'- (6.54) 



So this gives us: 

1 1 



xw/l 11c Sl,S2,Sj^l,S2/,Cl,C2,C^l,C2l 



2/,Cl,C2,C^l,C2l 



x^JiM£dS=5,(fc)5,(fc')5M'^2,2', (6.55) 
^x(l - x)x'{l - x') 



where: 



= IGTrV^^'^Pi - Pi')<^ci,c,,5.i,.„- (6.56) 

Our choice of basis functions guarantees finite kinetic energy, so we can drop the 
infrared regulator, 9^. Next, we can use: 



Cl,C2,Cj/,C2/ Cl,Cj/ Cl 

and 



Sl,S2,Si/,S2/ 



E X^"^x:r-^'<5,,.^,<5,„,^,A*_,^_3^(0)A,,_,^,_,^,(0') 

- V^--.i-.2(0)^;'-.,-.2('/'')- (6.58) 

Finally, we will also use: 

(Pl - Pl')^^^^ (P2 - P2') - ^(P]^ - Pr)^^^^ (Pl± - Pl'±)<^^^^ (P2 - P2') 

= P+x(l - x)6{x - x')6^^\[x - x']V± + [k± - k'^])5^^\[pi + p2] - \pv +P2']) 
= V+x{l - x)5(x - x')5^''\k^ - k'j)5^^\V - V). (6.59) 
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This leaves us with: 

167r3p+5g,g,5(3)(P - V') J (fk^(fk'^dxdx'5{x - x')5^^\ki_ - k'^) 

= 167r3p+(5g,q/(5(=^)(^ - / d'^k^dx 

xB,{x)Biix)B,{k)Btik)A*_,^_J<P)Aj.^s,-sM 
= lQ7r'^V+6g^g'6jj>6^^\V -V) J dkkdxBi{x)Bi,{x)Bt{k)Bt,{k). (6.60) 

There are three parameters with dimensions of mass in the meson problem, the par- 
ticle mass, m, the relative transverse momentum, k, and the cutoff A. We want to 
write the transverse momentum in units of the cutoff so the transverse B-splines will 
have a dimensionless argument. For the overlap in Eq. ( |6.6CI| ) to have the correct 
dimensions, the B-splines of the dimensionful parameter, k, must have dimensions 
. If we introduce the dimensionless B-spline Bt{k) with dimensionless argument 
k so that: 

^B^ik) = B,{k), (6.61) 

we can let 

dkkBt{k)Bt'{k) ^ dkkBt{k)Bt'{k). (6.62) 

This change adds extra notation that complicates expressions, although the expression 
has not changed. We drop the 'hat' from Bt and k and all B-splines are dimensionless 
and all transverse momenta are in units of the cutoff. We now write the overlap: 

{q',l',t',f\q,l,t,j) = 

lQ7r^V^6g,q>6jj'6^^\V -V) J dkkdxBi{x)Bi,{x)Bt{k)Bt'{k). (6.63) 
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This also shows that the B{k) have dimensions of inverse transverse momentum. 
When we make the change of variables: 

2 



k 



y + 1 



(6.64) 



we find: 



{q',l',t',j'\q,l,t,j) 



dy 



-iil + y)^\l+y J 
xBt{y)Bt,{y) f dxBi{x)Bv{x). (6.65) 



Finally, it is conventional to drop the factor IQn^V^ S'^^^ (V — V) as it appears in all 
matrix elements. So this gives us the final form: 



^ 1 j Bt{y)Bt^{y) £ dxBi{x)Bi,{x). (6.66) 



6.7 qq Eigenvector Normalization 



For our state to have the normalization in Eq. ( |6.17| ), it is necessary for the 
momentum part of the wave function to obey the relation: 



d^k^dxe{x)e{i -x)Y, WJls^i^^ 



1, 



(6.67) 



where 



4 ki — l kt CO 

7^ E E EE 

Xi?^,L^^,T:'a'AK^)A(^)^*(^)^HA:)^:'(0)^a(0),(6.68) 



and 



(6.69) 



which is derived in [27] based on symmetry about the 3-axis (there is a subtle differ- 
ence for quarks, since the quarks have spin ±| instead of ±1). After some algebra, 
we find: 



^ = <! E E ^J^:0^' \ , (6.70) 

q=l l,t,l',t' 



where 



Ofi, = / dxBi{x)Bv{x) / dkkBt{k)Bti{k), (6.71) 
Jq Jo 

and the i^^JJ are the eigenvector elements determined by diagonalizing the Hamiltonian 
matrix. The calculation of also leads to the expression for the spin-averaged 
dimensionless probability density: 



Uix,k/A) = 2nAkJ2\KJ^,k, 

At. 4 ^ ^ 

= ^ E E E R^ii!;i>;B,{x)B,{x)B,{k)B,{k) 

q=l l,t I't' 

2 

^ Y.KtBi{x)Bt{k) . 

Lt 



(6.72) 

6.8 Testing the Basis: Glueball Masses 



The calculation of the glueball spectrum in |2g, ^ uses oscillating, orthogonal 
basis functions that may be a source of round-off error in the numerical calculation 
of matrix elements. In addition it is likely their oscillations slow the convergence of 
the Monte-Carlo integration. Our approximate meson states use B-spline basis func- 
tions (Chapter |^) for the longitudinal- and transverse-momentum degrees of freedom 
because we believe they will avoid round-off error and their simple structure may 
speed the convergence of the Monte-Carlo integration. Before we commit to deriving 
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the necessary meson matrix elements in our basis, we check to see if there is any 
improvement in the calculation of the glueball spectrum using our longitudinal- and 
transverse- momentum basis functions defined in Sections |6.4. 1| and |6.4.2| , respectively. 

We briefly describe the calculation of the glueball mass spectrum in pure-glue 
QCD. Then in Section |6.9|, we write the orthogonal functions used in the previous 



glueball calculation, followed by a comparison of our results with those in [^, [27 



The second-order pure-glue QCD Hamiltonian matrix with gg external states is 
similar to the second-order full QCD Hamiltonian matrix with qq external states. 



They are so similar that the same diagrams and interactions shown in Section ^ 
with quarks and antiquarks replaced with gluons, describe the glueball system. 

The coupling is the only fundamental parameter in pure-glue QCD. Thus we only 
need to determine the coupling and cutoff that most accurately reproduce the glueball 
spectrum.]^ We make several calculations with different values of the coupling. Each 
calculation produces eigenvalues of the renormalized invariant-mass operator, A^^(A). 
These eigenvalues are used to determine the coupling and cutoff. Our procedure 
should produce approximate cutoff-independent results. The optimal value of the 
coupling is determined by first finding the range of couplings where the results are the 
most cutoff-independent and then within this range, finding the value of the coupling 
that reproduces the expected degeneracies for different values of j, the projection of 

the total spin along the 3-axis. 

-•^^There is no real experimental data for glueballs. Lattice gauge calculations are the most trusted 
source of data for the glueball system, so we compare our results with Ref. |40[. 
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Once the coupling is determined, we find the cutoff by fixing the first excited 
glueball state mass using the relation: 



2 

\'2 ''measured fa 

where "^measured the measured mass-squared for the state. 
6.9 Glueball Results 

In this section we show that our results for the glueball mass spectrum is consistent 



with [2(:, 27|. We begin by simply stating the basis functions used in the previous 



calculation (Section |6.9.1|) followed by a presentation of the results. We use kl to 



designate the B-spline parameter k for the longitudinal B-splines and k2 for the 



transverse B-splines. Section |6.9.2| shows that all desired results have converged when 

the basis parameters are kl = k2 = 6. This means there are a total of 7 transverse 

states, 8 longitudinal states (4 symmetric and 4 antisymmetric) and 4 spin states, for 

a total of 112 basis states. All plots except for the ones showing this convergence use 

these 112 basis states. Some matrix elements are zero by definition due to spin- flip 

restrictions.0 The four spin states split the Hamiltonian into sixteen blocks for each 

combination of incoming and outgoing spin. If j, the total spin along the z axis is 

not zero, there are only two blocks that are zero. However, if j is zero, then there are 

eight more blocks that are either zero or identical to another block. Thus, only six of 

the sixteen blocks may need to be determined.^ 

We want to make clear that although we do not include as many results as are 

given in we are able to reproduce all of the previous results. 

^°For example, there are no up-up to down-down interactions. 

^^The Hamiltonian is also Hermitian which further decreases the number of unique elements. 
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6.9.1 Orthogonal Basis Functions 

The difference between the basis functions we use and those found in [EO, BTIl are 



the functions used to represent the longitudinal- and transverse-momentum degrees 



of freedom. We use functions based on B-splines that are defined in Sections |6.4.1 



and |6.4.2| . The orthogonal longitudinal functions used by Allen are given by: 

I 



mx) = [x{l - x)r E (6-74) 



m=0 



where the A; ^ are determined by Gramm-Schmidt orthogonalization. The orthogonal 
transverse functions are functions of the relative transverse momentum k with the 



range [0, oo). They also depend on the parameter d with dimension 
to help minimize the ground state energy. They are given by: 



which is fixed 



t 

Tt^'^\k) = rfexp(-k2d2)Et^t,skM^ (6.75) 

s=0 

ft{kd) = -M''\k). (6.76) 
a 

6.9.2 Convergence Testing 

Section |3.1| introduces the ideas of convergence used in this section. The first 
important point is that the real ground state energy is always lower than the lowest 
eigenvalue of the approximate Hamiltonian. Conversely, the lowest eigenvalue pro- 
vides an upper limit to the real ground state energy. The second important point 
is that for a "standard" set of basis functions, as more functions are added, a bet- 
ter approximation to the real state is achieved, lowering the eigenvalue at each step. 
However, each time a new B-spline is added, the other splines in the set are rear- 
ranged. This makes it possible for eigenvalues to increase. This is typically a small 
effect, and has not increased the eigenvalue for more than one function addition. 
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This section shows how the eigenvalues of the IMO converge as we increase the 
number of transverse and longitudinal basis states. In the plots, the number of 
transverse basis states is A;2 + m + 1 — 3 where three are removed to ensure finite 
kinetic energy (see Section 6.4.2| ). The total number of longitudinal basis (symmetric 



plus antisymmetric) states is m + A;l + 1 — 2 where two that do not have the proper 
a; — > 0, 1 behavior, are removed. Increasing the number of longitudinal states has 
very little effect on the eigenvalues, and increasing the number of transverse states 
leads to rapid convergence for low-lying states. This suggests the B-splines in the 
longitudinal direction are very efficient basis functions for that direction. Figures |6.1| 
through |6.3| show how the eigenvalues of the mass-squared operator converge as kl 
and k2 increase when a = | for the ground state through the second excited state, 
respectively. These figures use m = 3, a = | and j = 0. This choice of parameters is 
arbitrary because we just want to see how many B-splines are needed in the transverse 
and longitudinal directions to give converged results. The "2% error" point plots the 
lowest eigenvalue and has 2% error bars. Thus any lines within these error bars are 
consistent with the lowest eigenvalue. 

6.9.3 Determining the Cutoff 

In this theory, the only dimensionful parameter is the cutoff, A. The cutoff can 
be determined by fitting the mass of one state. Explicitly, if we fix the cutoff using 
the i*^ state we find: 

2 

A 2 ''measured 

^ = (jmrn^ ^ ^ ^ 

^^The eigenvalues are calculated to only 2% accuracy. Any combination of kl and k2 that yield 
results within 2% of the lowest eigenvalue produce converged results. 
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n=0 




Figure 6.1: Ground state eigenvalues of the IMO for different numbers of longitudinal 
and transverse basis states with zero spin and the coupling equal to one-half, kl 
and k2 are the B-spline parameter k for the longitudinal and transverse B-splines, 
respectively. 

n=1 




0.59 I 

23456789 

k2 



Figure 6.2: First excited state eigenvalues of the IMO for different numbers of longi- 
tudinal and transverse basis states with zero spin and the coupling equal to one-half. 
kl and k2 are the B-spline parameter k for the longitudinal and transverse B-splines, 
respectively. 
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n=2 




Figure 6.3: Second excited state eigenvalues of the IMO for different numbers of 
longitudinal and transverse basis states with zero spin and the coupling equal to one- 
half, kl and k2 are the B- spline parameter k for the longitudinal and transverse 
B-splines, respectively. 



where ?Timeasured the measured mass-squared for the state. There is no real exper- 
imental data for glueballs, so we use lattice results [EO] for the "measured" values. 



Figure |6.4| shows how the coupling and cutoff are related when we fix the first excited 
state. 



Figure |6.5| shows that when we fix the first excited state, for a = .5, fcl = 6 and 
k2 = 6 the value of the cutoff has converged. 

6.9.4 Masses versus Coupling 

In this section we show how the masses depend on the coupling, a, when the 
first excited state is fixed. However in these plots rather than showing the mass's 
dependence on the coupling, we show its dependence on the cutoff. These plots help 
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a vs. A 

1 I f 



= 0.5 




1 1.2 1.4 1.6 1.8 2 2.2 2.4 

A/m 



Figure 6.4: The coupling a as a function of the cutoff A, in units of the first excited 
state mass when the first excited state is fixed. 




Figure 6.5: The cutoff for a few values of k2 when a = .5. This shows the result is 
converged. 
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Figure 6.6: Masses for j=0 states as a function of the cutoff when the first excited 
state is fixed. 



determine for what range of cutoffs the masses have a small cutoff dependence to 
determine where the theory is least cutoff dependent. In agreement with Allen's 



results ||2^, |2^, we find the slowest cutoff dependence for .5 < a < .7. 

In figures |6.6| through |6.8| the masses are plotted as a function of the cutoff when 
the first excited state is fixed. This is equivalent to plotting the masses as a function 
of the coupling. 

6.9.5 Glueball Spectrum 



For completeness we present the glueball spectrum from [^] in this section. Our 
results are identical to within statistical errors to these previous results. 

States are identified by J^''^ where J is the total spin of the state, j is the spin 
projection along the 3-axis, P is the parity of the state, and C is the state's charge 
conjugation eigenvalue. An asterisk in the state notation next to the value of C 
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Figure 6.7: Masses for j=l states as a function of the cutoff when the first excited 
state is fixed. 











j=2: masses vs. A 
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Figure 6.8: Masses for j=2 states as a function of the cutoff when the first excited 
state is fixed. 
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State 


M/Mo++ 


Lattice 


0-+ 


1.38 ±0.02 


1.34 ±0.18 


2++ 


1.58 ±0.01 
1.58 ±0.02 
1.11 ±0.01 


1.42 ±0.06 


2++* 


1.70 ±0.01 
1.68 ±0.02 
1.62 ±0.02 


1.85 ±0.20 


0++* 


1.77 ±0.02 


1.78 ±0.12 



Table 6.1: The glueball masses from p^, which are identical to our masses within 



statistical errors, compared to an average of lattice results from a number of different 
calculations |^ . Masses are in units of the mass of the 0"*"^ state. The uncertainties 



are only the statistical uncertainties associated with the Monte Carlo evaluation of the 
matrix elements of A^^(A). The three values of the masses for the 2"*"+ and 2++* states 
for our calculation correspond to j = 0,1,2. We use a = 0.5, with 8 longitudinal 
(4 symmetric and 4 antisymmetric) basis functions, 7 transverse-magnitude basis 
functions, and 4 spin basis functions, for a total of 112 basis functions. 



denotes an excited state with the given quantum numbers. Finally, we need to dis- 
tinguish states with identical J's and P's and different j's because we do not have 
manifest rotational symmetry. If J = 0, we omit the subscript j in the state notation. 



Table |6.1| lists the glueball masses for a = .5, in units of the mass of the ground 
state (the O"*""*" state). It also gives the average of lattice results from a number of 



different calculations for the sake of comparison The uncertainties listed are only 
statistical errors due to the Monte-Carlo integration. Finally, the three masses listed 
for the J = 2 states correspond to j = 0, 1, 2, where j is the spin projection along the 
3- axis. 



99 



Figure 6.9: Wavefunction for the 0"''+ glueball with a = |. 

6.9.6 Wavefunction Plots 

In this section we plot the spin-averaged probability density for the glueballs which 
is defined in Eq. ( |6.72| ) for states some of the low- lying states with j = and j = 1. 



We only show a few functions since they are similar to those found in ||2^, . The 
variable is the longitudinal momentum fraction for one of the gluons, and 

the variable is the relative transverse momentum of the state in units of the 

cutoff. 



Figures through |6.12| show the probability density for four of the the low-lying 
glueball states. Each wavefunction is built from 8 longitudinal (4 symmetric and 4 
antisymmetric) basis functions and 7 transverse functions. Including the four spin 
states, there are 112 basis states. 
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Figure 6.12: Wavefunction for the 2q ^ glueball with ol = \. 
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CHAPTER 7 



Matrix Elements of the Approximate Meson States 



We calculate all of the matrix elements of the approximate IMO (introduced in 
Section p.6|) necessary for a second-order calculation. The matrix elements are divided 
into four types of contributions; kinetic energy (KE), self-energy (SE), exchange (EX) 
and instantaneous exchange (IN). The four contributions are combined to explicitly 
cancel divergences (discussed in Section ^.7] ). With this division the full expression 
can be written: 



{q',l',t',f\M'{A)\qJ,t,j) = {q'J',t',j'\M'{A) 

+ {q',l',t',f\M\A) 

+ {q',l',t',f\M\A) 

+ {q',l',t',f\M\A) 

+ {q',l',t\f\M\A) 



qJ,t,j)KE 

9) ^) ^) j)ex 



IN+EX 



(7.1) 



where the free-state functions \q,l,t,j) are defined in Eq. ( |6.47| ). In the remainder 
of this Chapter we determine the matrix elements that need to be calculated by 
numerical integration. The integrals are divided into two groups: two-dimensional 
and five-dimensional. 
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7.1 Two-Dimensional Integrals 
7.1.1 The Kinetic Energy 

The two-dimensional integral for the kinetic energy is very similar to the overlap 
integral. The only difference in the derivation is: 



hq2\m2) 



x{l — x) 



(gi92l9l92)- 



(7.2) 



Thus, we have: 



\ 2 



dyk 



1 + 



1 + 



1 {l + vkf \l + yk 

1 , BAx)Bv{x 



1 



Bt{yk)B,{yk) f dx^^ 
Jo xii — 



X) 



(7.3) 



7.1.2 The Self Energy 



We need to separate Eq. ( ^.52|) into divergent and finite pieces: 



(M4|^sE(A)|gig2) 

<^l,3^2,4^2 



247r2 



z dz {l + z^ 



1 — z \ z 



^ [erf (7(x, z)) + erf (7(1 - x, z)) - 2] 



[7(x, z)Ei (1, 72(x, z)) + 7(1 - X, z)Ei (l, 7^(1 - x, z))] }] . (7.4) 



The divergent part of the self-energy is contained in the term: 

1 + z^ 



^1^3^24^2 

247r2 



J I — 7 (^2V27r 



-2V2^ ' ' A 
247r2 



21og(e) 



(7.5) 



So now we can split the self-energy contribution into a finite and divergent piece: 
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z dz ( 1 + z"^ 



247r2 



1-Z 



27r [erf (7(x, z)) + erf (7(1 — x, z))] 



+V8[j{x,z)Ei{l,-f\x, z) 

+7(1 - X, z)Ei (1, 7^(1 - X, z))] } + 3v^] , (7.6) 



{q3q,\vi\A)\qiq2) 



D 



247r2 



A^4V27rlog(e), 



(7.7) 



with 



{q'J\t',f\M\A)\qJ,t,jr^^ = gie 



2 -k-^ 1^1,1 I If +1 ■/|T/(2) 



(g',r,t',/|\/sW(A)|g,/,t,jr, (7.8) 



(g',r,t',/|A^^(A)|g,/,t,j)^E = ^7le 



(g',/',t',/|V^SE(A)|g,/,t,jr- (7.9) 



We find for the finite part (see Eq. |6.47 ): 



1 



X 



Mi(x)5Kx)I(x) / dqqBM)BAq). (7.10) 



where 



I(x) = 3v^+ / ^|l±^v^[erf(7(x,z))+erf(7(l-x,z))] 
J 1 — z [ z 

+V8 [7(x,2)Ei(l,72(x, z)) 

+7(1 - X, z)Ei (1, 7^(1 - X, z))] } . (7.11) 
The divergent part of the self-energy is identical, except: 

I(x)^4v^log(e), (7.12) 

giving 



(g', t\ 3\M\A)\q, l,t,j)^^ = (27r) i log(6) glV^S^'^ {V - V) 

/■l poo _ 

X5g,,'5jj, dxBiix)Bi,{x) dqqBt{q)Bt>{q). (7.13) 
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7.2 Five-Dimensional Integrals 
7.2.1 The Instantaneous Interaction 

The matrix element of the instantaneous interaction in the approximate basis is 
given by: 



{q',l',t',f\M\A)\q,l,t,j)i^ = gle-^-'^'^' X 

1 1 ^-^ /■ d?qdx f dPpdy 



E 



S,(?/)SKa:)5,(g)S,,(p)A*,_,3_,^(0O^,-.,-..(0)(g3g4|Vig)(A)|gig2), (7.14) 



where the matrix element of the reduced interaction is defined in Eq. ( |5.44|) . This 
can be reduced to: 



(g',r,t',/|A^^(A)|g,/,t,j)iN 
32 



g'\P~^6{V — V)6jj' J d'y^dqqdx J dppdyfteO, 
xd{x -y- e)Bv{y)Bi{x)Bt{q)Bt,{p) 



e ^Fi 



{x - yY 

X COs([j - 1]7<^) + 5g,2V,2 COs([j + 1]7<^) 

+ ((^9,3^,3 + '59,4V,4) COs(j70)] . (7.15) 



Finally, as discussed in Section |2.7|, we split the instantaneous interaction into two 



parts, a part above the cutoff that vanishes in the limit A ^ oo, and a part below 
the cutoff that remains in this limit: 



{q\l\t\j'\M\k)\q,l,t,j)t^ 



(g',/',t',/|A<2(A)|g,/,t,j)fN 



1 — e^^^ ^FK^KI 



x{q'J',t',f\M\A)\qJ,t,jU, (7.f6) 



e ^FK^KI 



x{q'J\t',f\M\A)\qJ,t,jU. (7.17) 
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7.2.2 The Finite Part of the Exchange Interaction 



The matrix element of the finite part of the exchange interaction is: 



j',l',t',f\M\A)\q,l,t,j)Ex = gle-''''^'-' X 
1 1 ^ f (Pqdx f (fpdy 



E 



-R R ,S S vSlS2 S3S4 y 



Bi\y)Bi{x)Bt{q)Mp)A}.s,-sM')^J-s.-^^^^^ (7-1^ 



where the matrix element is given in Eq. ( p.41| ). After many pages of algebra the 
contribution from the finite part of the exchange interaction is given by: 



(g',/',t',/|A^^(A)|g,/,t,j) 



EX 



d'j^dqqdx J dppdy 



1 Biiy)Bi{x)Bt{q)Bt>{p) 
x-y x{l-x)y{l-y) 



where the Sq^q' are given by: 



X 



a(1) a(1) 



1 2A-4A<'' a'^' 



^j,j'^q,q'^ (7-19) 



2,2 



S: 



3,3 



4,4 



X-y 



q^yil -y){l- 2x) - p'^xil - x){l - 2y) cos([j - 1]70) 



1 



x-y 



-pq 



y\l - 2x) -x^l- 2y) - 2y{l - x) cos(j70) 



qS{l -y){l- 2x) - p'xil -x){l- 2y) cos([j + 1]7^) 



-pq 



y\l - 2x) -x'il- 2y) - 2y{l - x) cos(j70) 



- cos(j7<^) [pgcos7<^(l - x - y + 2xy) 
1 — X — y 



x-y 

cos(j7<^) [pg cos 7,^(1 - x - y + 2xy) 
1 — X — y 



— y)q^ — x{l — x)p^^ — m^{x — y)' 



x-y 



y{l — y)q — x{l — x)p ) + m (x-y) 
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'S'l,3 = 

Ss,! = 

'S'1,4 = 

'S'4,1 = 

'S'2,3 = 

'S'3,2 = 

'S'2,4 = 

'S'4,2 = 

'S'3,4 = 



[(1 - 2|/)gcos(j7^) - (1 - x - |/)pcos([7' - 1]7^)] , 



im 



[(1 - a; - i/)gcos([j - Vfi^) - (1 - 2a;)pcos(j7</,)] , 



im 

72 
im 

im 
72 

72 

im r 

71 

im 



-(1 - 2y + 2y^)gcos(j7^) + (1 - y - x + 2xy)pcos([j - 1]7^) 
(1 - 2a; + 2a;^)pcos(j7<^) - (1 - x - y + 2a;y)gcos([j - 1]70) 
[(1 - 2|/)gcos(j7</,) - (1 - a; - |/)pcos([7' + 1]7<^)] , 
[(1 - a; - y)gcos([j + 1)70) - (1 - 2x)pcos{j-i^)\ , 
(1 - 2y + 2y^)5cos(j70) - (1 - y - a; + 2a:y)pcos([j + 1)70) 
(1 - a; - y + 2a;y)gcos([j + 1]'^^) - (1 - 2a; + 2x'^)pcos{j'y^) 



V2 

V2pq sin 7^ sin(j7<^) (1 - a; - y) , 



'S'4,3 = V2pq sin 7^ sin ( j 7,^) ( 1 - a; - y ) . 
2.3 The Divergent Part of the Exchange Interaction 

The complete contribution from the divergent part of the exchange diagrams is: 
{q',l',t',f\M'{A)\q,l,t,Mx- 

X J d-f^dqqdx J dppdy^-^^^^^^^Bi,{y)Bi{x)Bt{q)Bt>{p) 



x-y 



X 



a(i) a(i) 



1-e 



p^ + q^ - 2pgcos(7^) \ 
{x - y)2 J 



X C0S([7 - 1]7<^) + 5g,2(^q',2 C0s([j + l]-^^) 

+5q,z5q',Z C0s(j7<^) + (5q,4(^q',4 COS(j7^)] . (7.20) 
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7.2.4 Combining the Divergent Part of the Exchange Interaction and the 
Instantaneous Interaction Above the Cutoff 

If we consider the contribution from the instantaneous interaction above the cutoff 
(Eq. [7.1 6| ) and the divergent part of the exchange interaction (Eq. [7.2(]| ), we see that 
in the limit x = y the two divergent contributions cancel. Thus, we can combine the 
two interactions and write them as: 



where 



(g',/',t',/|^^(A)|g,/,t,j) 



IN+EX 



j d'j^dqqdx J dppdy9^9^'9(x — y — e)e ^ ^^^^ 



xBi,{y)Bi{x)Btiq)Btip) 
1 



{x - yf 



1 2A-*A*^' a'^' 



1 \ (p^ + q^ — 2p- q 



.(1) 



x-y 



+ 2 



, (7.21) 



Sq,lSq>^l COs([j - 1]70) + Sg^2Sq',2 COs([j + l]-f^) 

+ ((^5,3^,3 + ^qA^Q'A) COs(j7<^). (7.22) 



In the limit \x — y\ the term in square brackets (inside the braces) vanishes, 
leaving the entire term finite in this limit so we no longer need the infrared cutoff for 
this term. 

7.2.5 Combining the Divergent Part of the Self-Energy and the Instan- 
taneous Interaction Below the Cutoff 

The combination of the instantaneous interaction below the cutoff and the diver- 
gent part of the self-energy is algebraically complex. The details of the calculation 
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can be found in Appendix O. We find: 



{q',l',t',j'\M\A)\q,l,t,j)^/ = -gl-S,,,V^6{V - V')9{x - y)e{x - y - e) 



where, 



1— e px—e 

dx 

X 



dy\og{x-y) / d'qd^pBi{x)Bt{q)Btip)Wq^^>e 



My) , B'Ay) B,{y) 



{x - yy 



x-y 



x-y 



(7.23) 



^KI 



rn?{x — y)"^ + [xp — yq)"^ 

xy{x - y) 
rn?{x — y) — yq^ + xp^ {p — q)^ 

I ' 



xy 



x-y 



A 



FK 



vr?{x — y)^ + [(1 — x)p — (1 — y)q\ 

(1 -y){x-y) 
rv?{x — + (1 — ?/)g^ — (1 — (p — qf 



[\-x){\-y) 



x-y 



A' 



KI 



a; 



[p — qf' + 
(x - y)2 ?/2 

+ p^ {p — q)"^ 



FK 



{l-yf {x-yY' 

However, the integrand has an apparent divergence as x — » y. The integrals over 
transverse momenta, however, integrate to zero in this case, kilhng the divergence. 
We can subtract a term that integrates to zero and that exphcitly cancels the false 
divergence from the above contribution. We make the following change of variables 
to aid in the calculation of the term we want to subtract: 



r = \{p + q), w=]^-^, (7.24) 

where 

r = r [cos ax + sin ay] , w = w [cos 5x + sin (7.25) 
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and 



(3 = a — 5, r ■ w = rw cos (3. 



(7.26) 



Then 



q = r + rjw, p = r — rjw, 



(7.27) 



and 



g = -i/r^ + r^w^ + 2™,/// cos/5, p = J + rjw'^ — 2rw^/r] cos (3. (7.2 



For convenience, we define: 



r± = Jr'^ + rjw'^ ± 2rwy/r]cos (3. 



Finally, 



q - p = r'^ — r]w'^, 



cos 7,^ 



2 2 

r — r^w 
r+r_ 



r 1 r_ 



(7.29) 



—2rWy/ri 
sin7(^ = — ^ — sinp, (7.30) 



are useful relations. The details of this subtraction are in Appendix 0. After making 
the subtraction, we have: 

{q',l',t',f\M\A)\q,l,t,j)i/ = 

-gl—Sjj'V^SiV-V) dxBi{x) r dyr]\ogr] f Sfd^w 
3 ^0 -^0 -J 



Bv{y) B[,{y) 



V 



V 



5n al 



-Bt{r)Btir)Bi,{y)^e-''^ [l-QAw^\j. (7.31) 

The integral is explicitly finite, but still sharply peaked around x ~ We can speed 
convergence of the numerical integration by smoothing the integral using the change 
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of variables: 



T] — xe 



dy = Tjdp, 



(7.32) 



and 



1 + w 



- 1, 



1 + r 



-1, 



1 +p 



- 1. 



(7.33) 



Finally, we get: 



■\B,F 



{q'J\t\f\M\A)\qJ,t,j)^^ 
128 



X 



dyr 



dvu 



X'qlogT] r wBi{x) 

X {e-(^^-+^?-%,,,,St(r+)St,(r_) 



-Mr)Mr)Biiy)^-^e-''-' 



1 - 64w' 



, (7.34) 



where, 



Afk = 



A' 



FK 



m^r) — yr^ + xr^ ^ ^ ^ 
xy 

m^r) + (1 - ?/)r^ - (1 - x)r^ 
(l-^)(l-y) 



V y 

w? + 4^^ 
(l-y)2 
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CHAPTER 8 



Numerical Issues and Parallelization 



With a problem of this nature the analytic calculation is only part of the solution. 
To get results that can be used to check the theory or make predictions, it is necessary 
to further refine the calculation for efficient numerical calculations. Furthermore we 
can obtain results more rapidly from a well-designed algorithm that makes the best 
use of computer resources. 

In this chapter we discuss some of the refinements that help make the Monte-Carlo 
integration more efficient. We also describe the algorithm and how it is parallelized. 

8.1 Numerical Calculation of Meson Mass Spectra 

In order to aid the numerical calculations, we will use the change of variables given 
in Eq. (^^) through Eq. (fTSOl) and Eq. (^^) through Eq. (^^). The complete 
change of variables gives: 



dy dq dp ^ 32?7^ / — 
Jo Jo Jo J -I 1 



dys dy^ dyr ,^ 



+ ys?J-i {l + y^fJ-i [l + yrf 
We divide out the common factor IQn'^V^S^iV — V) from all the matrix elements 
and use the following definitions for all of the integrals: 
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Wq,q' = ^q,l^q',l COs([j - l]-f^) + Sg^2Sq',2 COs([j + l]7</,) 

+ ((^9,3^,3 + SqA^q'A) COs(j70) , (8.2) 



A 



_|_ ^2 ^2 _|_ ^2 



FI 



y{l-y) x{l-x) 

xil-x)yil-y) K[^(l-^)-l/(l-y)] 



+77(1 — X — ?/) (^r^ + 



-2™ 1/7/ cos /3 [x{l — x) + y{l — y)]} 



m^ri - yr\ + xr'i ^ 

Aki = h 4w 

xy 

rjm'^ + T] (r^ + riw"^) — 2rw^ cos p(x + y) , 

= h Aw , 

xy 

m?ri + (1 — y)r\ — (1 — x)r'^_ 



^FK 



{l-x){l-y) 



rjw? + (r^ + r^w^) + 2rw^cos f3 {1 — x + 1 — y) , 

4w , 



{l-x){l-y) 



9 ' 

^2 



= (8.3) 



where we have rewritten some of the definitions to reduce round-off error. When 
calculating the cosine and sine functions of 7,^, we will need to use recursion relations 
since we can really only calculate sin70 and COS70 [see Eq. ( |7.30|) ]. 

8.1.1 The Finite Part of the Exchange Interaction 

From Eq. (|7l9| ) we eret: 
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{q\l\t'j'\M'{A)\q,l,tJ)l^^ 
1 



j^A / d-f^dqqdx j dppdy 



-A2 

e Fi 



Gtt^ " J ' J X — y 

Bi.{y)Bi{x)Bt{q)Bt,{p) ( 1 



X- 



.2 A I 



16 



37r3 



2a /^^ /^"^ dys dyyj 

1 



dyr 



^^^ Bi,{y)Bi{x)B,{yr^)BAyr.) ^^^-Al.r 1 



x{l -x)y{l-y) 



X 



1-e' 



(8.4) 



where 



and: 



yr± 



1 + r± 



1, 



(8.5) 



^1, 



5. 



r+i/(l - |/)(1 - 2x) - rlx{l - x){l - 2y) cos([j - 1]7<^) 



-T] r-r+ 



y\l - 2x) - x^l - 2y) - 2y{l - x) cos(j7^), 



2,2 



'S'3,3 — 



'S'4,4 — 



rly{l - 2/)(l - 2x) - r^_x{l - x){l - 2y)\ cos([j + 1]7^) 
-T) r^r+ [y\l - 2x) - x'^{l - 2y) - 2y{l - x)] cos(j7^), 
■ cos(j7<^) [r}r_r+ cos 7,^(1 -x -y + 2xy) 

— (1 — X — y) [y{l — y)r^ — x{l — x)r^^ — m^rf 
cos(j7<^) [77r_r+ cos 7,^(1 - x - y + 2xy) 

— (1 — X — y) {^{1 — y)r1_ — x{l — x)r'^^ + rri^r)'^ 
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>S'l,2 
5'2,1 
Sl,3 

5'l,4 
5'4,1 

'S'2,3 
'S'3,2 
>S'2,4 



4,2 



'S'3,4 
5'4,3 



0, 
0, 

imr] 

7f 

imt] 

'V2 

imt] 



imr] 

7f 



imr] 
imr] 
imr] 



imr] 



[(1 - 2y)r+ cos(j7<^) - (1 - x - y)r_ cos([j - 1]7<^)] 
[(1 - a; - y)r+ cos([j - - (1 - 2x)r^ cos(j7<^)] 
-(1 - 2?/ + 2?/^)r+ cos(j70) 

+ (1 - y - a; + 2a;?/)r_ cos([j - l]-^^)] , 
'1 — 2a; + 2x^)r_ cos(j70) 

-(1 - a; - ?/ + 2x?/)r+ cos([j - 1]7^)] , 
[(1 - 2y)r+ cos(j7^) - (1 - a; - y)r^ cos([j + 1]7^)] , 
[(1 - X - ?/)r+ cos([j + 1]7^) - (1 - 2a;)r_ cos(j7<^)] 
{l-2y + 2y^)r+ cos(j7^) 

-(1 - ?/ - X + 2x?/)r_ cos([j + 1]70)] , 
[(1 - X - ?/ + 2xy)r+ cos([j + 1]7^) 



— (1 — 2x + 2x^)r_ cos(j7(^) 
r]V2r_r+ sm'j^sm{j'j^){l - x - y), 
r]V2r_r+ sm-f^sm{jj^){l - x - y). 



8.1.2 The Instantaneous and Exchange Interaction 



Eq. ( [mD gives: 



116 



(g',/',t',/|A^2(A)|g,/,t,j)iN+EX = P+^^ie"^"'^- x 



16 

— - / d^^dqqdx I dppdyeA'Bv{y)Bi{x)Bt{q)Bt'{p) 



X 



1 - e 



X 



'^fk'^ki 



,(1) 



_|_ g2 _ 2p . q 



-^,9AhfldxJ^ d^, 



xe-^'^'Bi,iy)Biix)Btiyr^)Btiyr 



dys 



dVu 



dyr 



sw 



Aw' 



A 



1 , 

+ 2 



FK ^KI 



(8.6) 



8.1.3 The Instantaneous Interaction 



The final five-dimensional contribution comes from Eq. ( |7.23| ): 



(g',/',t',/|A^^(A)|g,/,t,j)f/ = -^gl^6,,.V^8{V -V) 



128 



1 |-27r rl 

dx I d'y^ 



JO 



dyp 



dyr 



1 {1 + ypyJ-i {l + yr^J^i (l + y. 



dyn 



7] logr/ r wBi{x) x 



+ B',{y) - 2Bv{y) {^fk^'fk + ^ki^'i 



KI) 



-Bt{r)B,{r)Bi,{y)^e 
V 



1 - QAw' 
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64 



2n 



dyp 



dyr 



dyu 



X log 77 r wBi{x) [e-^^'i-K+^'''K^Wq^q^Bt{r+)Bt'{r_) 



X 



B,{y) + r^B[,{y) - 2B,{y) (AfkA'^j, + AkiA'^j 



-Bt{r)Bt>{r)Biiy)6g^g,e-^'^' 
8.1.4 The Self-Energy 



1 - 64?^;' 



}■ 
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The self-energy is given from Eq. ( [7.10|) : 
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8.1.5 The Complex Hamiltonian 



One substantial difference between the meson and glueball Hamiltonians is that 
the meson Hamiltonian is complex. However, each matrix element is either real or 
pure imaginary, so there is still only one integral per matrix element. This means we 
can treat the problem the same way we treated the glueball numerically. Specifically, 
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ii q — q' or {q — 3 and q' = 4) or (g = 4 and q' — 3) the matrix element is real, 
otherwise it is imaginary. 

8.1.6 Reducing the Number of Matrix Elements to Calculate 

We can take advantage of certain features of the Hamiltonian to reduce the num- 
ber of matrix elements we actually need to calculate. First consider the exchange 
interaction since it contains the only contributions when q ^ q'. 

The exchange interaction contains non-zero contributions only when q ^ q'. The 
"main" part of the exchange diagram (everything except the ^'s) is obviously sym- 
metric when we switch final and initial states. This switch is accomplished simply by 
letting q q' , I ■(-^ I' , t ^ f, x y, q (r+) p (r_), and 7,^ — > —70- 

We see that the real S''s are also symmetric under interchange of final and initial 
states, and because they are real this part of the Hamiltonian is already Hermitian. 
The imaginary S^s are also antisymmetric under interchange of final and initial states. 
This guarantees the entire Hamiltonian is Hermitian and cuts almost in half the 
number of matrix elements that must be calculated. Of course Hermiticity is necessary 
to obtain real eigenvalues and it provides a check on the analytic work up to this point. 

More redundancies in the Hamiltonian can be found if we consider the case j — 0. 
First, consider the ^'s in this case: 



rly{l -y){l- 2x) - r'ix{l - a;)(l - 2y) cos{'y^) 



77 r_r+ y\l - 2x) - x''{l - 2y) - 2y{l - x) , 



'2,2 — ~ 



rly{l -y){l- 2x) - r^x{l - x){l - 2y) cos{j^) 



77 r_r+ y^l - 2x) - x'(l - 2y) - 2y{l - x) , 
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'S'3,3 = - [vr-r+ cos 7^(1 - x - y + 2xy) 

— (1 — X — y) (^y{l — y)r^ — x{l — x)r^^ — w?rf 
— [77r_r+ cos 7^(1 — x — y + 2xy) 

—(1 — X — y) {y{l — y)r]^ — x{l — x)r't^ + m^rf 

0, 



'J4,4 

'S'1,2 
'S'2,1 
'S'1,3 
5*3,1 
'S'1,4 
'S'4,1 
'S'2,3 
'S'3,2 
'S'2,4 
'S'4,2 
'S'3,4 
'S'4,3 



= 0, 



imr] 

71" 

imr] 

72 
imr] 

72 
imr] 

71" 

imr] 

71" 

imr] 

72 
imr] 

imr] 
V2 



[(1 - 2y)r+ - (1 - X - y)r_ cos(7<^)] , 

[(1 - a; - y)r+ cos{'y^) - (1 - 2a;)r_] , 
-{l-2y + 2y^)r+ + {1 - y - x + 2xy)r- cos{-f^) 
(1 — 2a; + 2a;^)r_ — {1 — x — y + 2xy)r+ cos(7<^) 

[(1 - 2y)r+ -{1-x- y)r_ cos(7<^)] , 

[{l-x- y)r+ cos(70) - (1 - 2x)r_] , 
(1 - 2y + 2y\+ - {1 - y - x + 2xy)r- cos(7<^) 
{1 — X — y + 2xy)r+ cos(70) — (1 — 2a; + 2a;^)r_ 



= 0, 



= 0. 



Prom this we see that Si^s — ^2,3, '^1,4 = —82,4, and §3^4 — 0. 

The only redundancy when q = q' is when q = l,q' = 1 and q = 2,q' = 2. So we 
have the foUowing symmetries: 

{2,l',t',0\M^\2,l,t,0) = {l,l',t',0\M^\l,l,t,0), 
(l,/',t',0|A^2|3,/,t,0) = (2,r,t',0|A<^|3,/,t,0), 
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{l,l',t',0\M^\4:,l,t,0) = -(2,/',t',0|M^|4,/,t,0), 
{3,l',t',0\M'^\4:,l,t,0) = 0. (8.9) 

8.1.7 Rotational Symmetry: j —j 

Reference shows that the basis \q,l,t,—j) is related to the basis \q,l,t,j) 
by swapping the states \l,l,t,j) and \2,l,t,j), and changing the sign of |4, 
Renaming the basis states and changing their phases has no effect on the eigenvalues 
of the matrix, so the eigenvalues of states with the same absolute value of j are the 
same. 

8.2 Hamiltonian Matrix with Parallel Processing 

Fundamental physics research is dominated by analytical calculations. However, 
as theories become more complex, comparing theory to experiment becomes more 
numerically intensive. Fortunately the computational performance to price ratio has 
been increasing at such a rate that many numerical problems previously considered 
impractical can be done on readily available desktop computers. 

Although computers have and will continue to get faster, algorithms written to 
run on parallel processors almost always perform better than their single processor 
counterparts. The amount of work that is required to write a program using such 
standards as the Message Passing Interface (MPI) may seem inhibiting, but with 
certain algorithms, like the one used to numerically determine Hamiltonian matrix 
elements, there can be a tremendous benefit. 
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8.2.1 Thread Independence 

Algorithms used in fields such as Lattice Gauge Theory use program threads that 
directly depend on each other. For instance, the lattice may be equally divided 
among the threads, and then boundary information is passed between them. The 
efficiency of this type of algorithm, whose threads communicate frequently, can be 
largely inffuenced by the network environment. However, when calculating matrix 
elements, each matrix element can be calculated almost independently of the others. 
Thus the communication between threads is very infrequent and does not greatly 
decrease the performance of the program. This fact allows our algorithm to run 
efficiently on many different machines over a large, inefficient network. 

The majority of our algorithm used to calculate the Hamiltonian matrix is em- 
barrasingly parallel, each worker process only being told which matrix element to 
calculate. However, there are parts of the program (Sec. 8.2.3| ) that use a parallelized 



version of the numerical integrator (Sec. ^.2.6|) to limit idle worker processes. 



8.2.2 Determining Accuracy of Matrix Elements 

The longer a program runs, the more important it is to make sure the source code 
is written well to reduce run times. Similarly, it is important that the program is not 
performing calculations that do not affect the result. 

Each matrix element in the Hamiltonian is determined by numerically calculat- 
ing a five- dimensional integral using VEGAS, an adaptive Monte-Carlo integrator. 
Monte-Carlo integrators are much more efficient then nested single dimensional in- 
tegrators for a large enough number of dimensions.^ However, the results converge 
^•^For three or more dimensions Monte-Carlo is generally faster. 
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slowly^. Thus, each matrix element should only be calculated as accurately as needed 
to produce results within the desired error. 

If all matrix elements are calculated to the same precision, the program will run 
for many times longer than necessary. Although it is dependent on the value of 
input parameters, to achieve two-percent accuracy in the eigenvalues only about ten 
percent of the matrix elements are important enough to require an error less than 
fifty percent. Unfortunately which matrix elements are important can not be known 
a priori. The method used to calculate each matrix element to within the required 



precision is discussed in Section 8.2.3 



8.2.3 The Simplified Algorithm 

The purpose of the algorithm is to calculate all Hamiltonian matrix elements to 
sufficient accuracy to give eigenvalues within a pre-determined error. Two questions 
must be answered before we can determine how accurately to calculate each matrix 
element. How should the error in the eigenvalues be determined from a matrix in 
which every matrix element has a known uncertainty? How does the uncertainty in 
a matrix element translate into errors in the eigenvalues? 

Before answering these questions we give a little bit of detail about how the 

matrix element is actually determined. Each matrix element has a contribution from 

a two-dimensional and a five- dimensional integral. The two-dimensional contributions 

are determined by nesting two one-dimensional integrations. This first contribution 

can be calculated easily and quickly to a much higher accuracy than the remaining 

contribution from the five-dimensional integral. VEGAS saves statistical information 
^''Monte-Carlo results converge like for N function calls. 
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about the integral, so if the desired accuracy is not reached, the routine can be re- 
entered, and the calculation will proceed from where it left off. If the routine did 
not have the ability to be restarted with previously obtained information, the desired 
precision would need to be known a priori. Since matrix diagonalization is a highly 
non-linear problem, understanding how each matrix element affects the answer is 
nearly impossible. 

With no reasonable analytic way of understanding how the uncertainties in the 
matrix elements translate to errors in the eigenvalues, we calculate the error statisti- 
cally. Each matrix element is independently and randomly varied about its average 
value using a gaussian distribution determined by the standard deviation, then a set 
of eigenvalues is determined. This is repeated to generate a large number of eigenvalue 
sets. Since we only expect the lowest several eigenvalues to be correctly determined 
by our renormalized Hamiltonian, we calculate the average and standard deviation of 
the lowest ten eigenvalues. The largest uncertainty in any of these lowest eigenvalues 
is defined as the error in the eigenvalues. Although this is a good statistical way to 
determine the error, it can be very cpu-intensive for large matrices, and thus is not an 
efficient way to determine the deviation produced in the eigenvalues for each matrix 
element individually. 

To approximate the deviation in the eigenvalues caused by a single matrix element 
we use a much less thorough method that only requires calculating the eigenvalues 
three times (versus a hundred or more). The eigenvalues are determined when the 
element under question is at its average value, and then when it is its average value 
plus and minus its standard deviation. The deviation is defined as the largest change 
in an eigenvalue when the element is shifted up or down by its standard deviation. 
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With this understanding of the deviation produced by a single matrix element, and 
the error in the eigenvalues from the entire matrix, the basic algorithm is straightfor- 
ward. After the two-dimensional integrals have been calculated, the algorithm cycles 
through each matrix element reducing the deviation produced by each element to 
a predetermined value. At the end of an iteration, the error in the eigenvalues is 
determined. If the error is less than the desired error, the program exits, otherwise, 
it cycles through the matrix elements again, this time requiring the deviation to be 
smaller than it was for the previous iteration. The program iterates through the ma- 



trix in this manner until the desired accuracy is reached. Figure shows the flow 
chart for the basic algorithm. 

8.2.4 Saving Data 

Saving basic data like the Hamiltonian and its error can be very useful because 
output can be rapidly regenerated or generated in a different format if needed. How- 
ever there is much more important data that can be saved. The most cpu-expensive 
data generated while running this program is the statistical information generated by 
VEGAS for each integral. 

In its basic form VEGAS only saves the information for the one integral it is 
calculating, and discards that information when it moves on to the next integral. 
The problem is that the same integral may need to be calculated more accurately 
in a later iteration, thus the integration would have to start from the beginning. 
Although saving the VEGAS information for each matrix element uses a lot of storage 
(30 integers, 1649 double precision), it is important to save all of it. If there is not 
enough memory available, the algorithm uses a scratch file to store the VEGAS 
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Simplified Algorithm for the 5D Contribution to a Matrix Element 
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Figure 8.1: Flowchart for the simple algorithm. The dashed line can be ignored to 
understand the basic flow, however in the parallelized version (Sec. ^.2.5| ) it separates 
the dispatcher's flow (left) from the worker's flow (right). 
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information for the entire matrix. During most of a new calculation 70 — 90%)^, 
there is little slowdown from disk access. 

Despite the performance gain in modern computers, they are still subject to 
crashes. Also, if the code is running in parallel on many different computers, the 
chance of the code stalling due to a single computer failure increases. This prob- 
lem can be remedied in a straightforward manner by occasionally saving the VEGAS 
information from memory to a file which can later be read in for code restart. 

One final benefit of saving data, in particular the data for the five-dimensional 
integrals, is that the coupling is an overall factor multiplying all five-dimensional 
integrals. If the parameters vary only in the coupling from a previous run, the old 
matrix element estimates can be used as a starting point, which speeds searches 
through our parameter space. 

8.2.5 Parallelizing The Algorithm 

We first define the nomenclature we use when discussing the parallelized algorithm. 
The Message Passing Interface (MPI) uses one dispatcher process and several worker 
processes.^ In this algorithm the dispatcher is the organizer. It sends jobs for the 
workers to do and receives their results. In addition the dispatcher process does not 
use a lot of computer cycles and therefore does not require its own processor. There is 
an important distinction between processes and processors. The processor(s) is(are) 
the computer's physical cpu(s), whereas the processes are simply jobs running on 

the machine. For example, on a particular machine, we typically run four processes 

^^A calculation that is restarted or one that uses saved data may go through several iterations 
reducing the deviation until there are integrals that need more precision. Until this point is reached, 
the VEGAS information is accessed frequently. If the information is stored on disk instead of 
memory, there is a significant but temporary slowdown in the calculation. 

^^There is nothing in MPI that forces the distinction between the dispatcher and the worker. 

127 



on two processors. One process is the dispatcher and the remaining three are all 
workers.^ 

There are many different ways this problem can be parallelized. For example, 
the matrix can be divided equally between the processors, each responsible for one 
submatrix. Another method is to let each process always work independently, moving 
from element to element until the desired error in the eigenvalues is reached. These 
two methods are examples of something that is too simple, and too complex, respec- 
tively. The first method assumes that each section of the matrix is equally important. 
However, the far off-diagonal elements are very small and calculating them to within 
more than an order of magnitude is often unnecessary. The second approach requires 
a very complicated method of tracking which process is working on which matrix 
element, trying to reach a particular deviation. Also it is not clear that this method 
would be efficient since only a handful of matrix elements are very important, and 
they take a lot of time to calculate to sufficient accuracy. 

The philosophy we use in parallelizing the program is one of efficient simplicity. 
The method is simple because the calculation is broken up into steps and because 
no process continues to the next step until all other processes are ready to proceed. 
Efficiency is achieved by having processes that have finished help the processes that 
are still working. 

Although the entire code is parallelized, we only discuss two sections in detail. The 
first is the parallelization of the error calculation and the second is the parallelization 
of the 5-d integral calculation. At the end of each iteration (reducing the deviation 

in the five-dimensional integrals) the error in the eigenvalues needs to be calculated. 

^^The reason we run more than one worker process per processor is to artificially slow the workers 
on faster machines so the cluster is more homogeneous. 
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However, particularly for large matrices, this calculation can take a long time because 
it requires the complete eigensystem analysis. If only the dispatcher is used to do 
the error calculation, it may take several minutes, leaving all of the worker processes 
idle. The error calculation is parallelized by having each worker process generate an 
equal number of eigenvalue sets, then sending the results to the dispatcher, who then 
determines the averages and standard deviations. 

The part of the calculation that benefits the most from parallel processing is the 
calculation of the five- dimensional integral contribution to each matrix element. Each 
matrix element can be calculated independently so there is almost no overhead due 
to communication between processes. A worker is given an assignment from the 
dispatcher. It then calculates it and returns the result, the process repeats with a 
new assignment, and an updated Hamiltonian. Figure the simple algorithm, also 
shows which part of the flow is executed by the dispatcher process (left of dotted 
line) and the worker processes (right of dotted line). The only time this is not true is 
at the end of each iteration. Once all of the matrix elements have been assigned, as 
individual worker processes finish, the dispatcher directs them to help another worker 
process complete its assignment. The implementation of this "helping" algorithm is 
the first important and difficult addition to the parallelized algorithm. 

If there is only one matrix element that requires a long time to calculate, then as 
the other workers finish, they sit idle waiting for the last element to be determined. 
This is a major drawback of the simple parallelization (independent matrix elements) 
since, except at the end of each iteration, the program can efficiently run using many 
processes. Thus, all but one process can be sitting idle while the last one continues to 
work. This can be remedied by having the dispatcher keep track of who is working and 
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who is available to help. This increases the communication between the dispatcher 
and the worker at the time of assignment. The worker must be notified if it will be 
helping another process, or if it will get help, and who will be helping it. When the 
worker is finished it must let the dispatcher know if it was helping or getting help, 
etc. Figure ^.2| shows the fiowchart for a worker process with the option to help other 



processes. 

The last complication that requires additional communication between processes 
is saving VEGAS information. Not only should the VEGAS information be saved for 
the completed matrix elements, it should also be saved for calculations in progress. 
Otherwise, information from an integral that has consumed hours of cpu time might be 
lost. Thus the worker may also send updated VEGAS information, without sending 
a result. It is non-trivial to work out the communication between the worker and 
dispatcher process to allow for a wide range of information transfer. 

8.2.6 Parallelized VEGAS 

The parallelization of the numerical integrator is done in an attempt to limit 
the total number of communications between the worker in charge of a particular 
integral, and the processes helping it. Also, we do not want to get into the details 
of the adaptive nature of the integrator but do want a simple way of speeding the 
convergence of the integral. 

The worker calls VEGAS with, among other things, the number of points to use 
and the number of iterations to attempt. After the desired number of iterations is 
completed, the routine statistically combines each of the iterations to get the result 
and its uncertainty. Thus, it is very straightforward to include information from the 
helper routines as extra iterations in the main worker. 
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Helping Algorithm For a Worker Process 




Figure 8.2: Flowchart for a worker process during the five- dimensional calculation, 
including the steps to help other processes. The double-lined arrows represent com- 
munication to and from the dispatcher. 
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The information that is lost using this method is that which helps refine the grid 
to make VEGAS an adaptive integrator. However, according to the documentation 
|1[ , the grid does not get significantly refined after a few iterations of a small number 



of points. Thus discarding the grid information should be unimportant. 
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CHAPTER 9 



Results 



In this Chapter we present our fully relativistic results for the meson mass spec- 



trum. The results are presented along with experimental data |^2[ and predictions 
from a light-front calculation using a similar renormalization method with a nonrel- 
ativistic reduction 

We first present important information about the experimental results we use. 
Our results begin with a brief description of the steps required to calculate the meson 
spectrum followed by two applications. We apply the method to the bb system as 
a check of the theory to verify that it agrees with nonrelativistic results and the cc 
system is investigated in detail. Finally, we show the spin-averaged meson probability 
densities defined in Eq. ( |6.72|) to show what type of wavefunctions are produced. 



9.1 Experimental Results 

Heavy quarkonia states can be produced by colliding electrons and positrons.0 
As the center-of-mass energy is changed, an increase in scattering events signals the 

existence of a resonance. The resonance is associated with a particle being created 

^^It is possible that some states couple very weakly to this production channel which could lead 
to missing states in the Particle Data Table. There is important ongoing work seeking new hadrons 
using new beams and targets. 
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and then decaying into other particles. The mass of the initial particle is determined 
simply by finding the center-of-mass energy at which the resonance occurs. Daughter 
particle masses are determined by subtracting the energy of emitted photons or pions 
from the parent particle. Identifying the particles generated in a detector is not 
simple. One must infer the quantum numbers by knowing the quantum numbers of 
the beam and the daughter particles that are actually detected. 

Although particle accelerators have produced a large amount of data that has been 
analyzed to give a plethora of information on spectroscopy of systems and transitions 
between particles, not all theoretically predicted particles have been seen. This fact 
is important when comparing new theoretical results to experimental data. 

The charmonium system contains two (IS) states, the J/ip and the 7]^. The rjc 
is not directly produced in colliding e^e~ experiments but it is observed through 
transitions from the J / if) and other heavier states. The bottomonium system is also 
expected to have two (IS) states but the 77^ has not been seen experimentally because 
the transition from the T is greatly suppressed relative to the J/ip •yrjc because of 
the different quark masses, charges and photon energies. 

9.2 Procedure 

It should first be noted that there is no simple way to guarantee correct state 
identification when performing this calculation because, for instance, the spin operator 
is dynamic on the light-front as is equal-time parity. To simplify state identification, 
the quantum numbers we use for each data set include the z-component of the total 

spin (j) and charge conjugation (C).0 The drawback is that there are fewer low-lying 

^^We can also (usually) determine the intrinsic spin of the system and the state's symmetry in 
the longitudinal momentum direction 
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states associated with each set of quantum numbers. For instance, the lowest five bb 
states contain two with negative charge conjugation {n = and n = 4), and three 
with positive charge conjugation (n = 1,2,3). We must use the C = + states to try 
to determine the mass ratios because we need at least two ratios to fix m and a.Q 

There are two fundamental QCD parameters that need to be fixed, and the scale 
(cutoff) needs to be determined. The calculation produces eigenvalues of the renor- 
malized invariant-mass operator, Ai'^{m^, A). These eigenvalues are used to deter- 
mine the fundamental QCD parameters and the cutoff. The steps to determine these 
values are described next. 

Eigenvalues are generated for a set of a and ^. If the renormalization was entirely 
non-perturbative and the Hamiltonian was calculated to all orders, the spectrum 
would be independent of the cutoff; however, we expect some cutoff dependence 
because we renormalize perturbatively and keep only two partons in our states. The 
renormalization is most trustworthy where there is little cutoff dependence. 

The first step in finding the correct value of a and ^ is to determine the cutoff 
by fixing one mass0 using the relation: 

2 

A 2 measured (QT^ 

~ {t\M^{A)\t)' ^ ' 



where '^measured taken from and i refers to the state used to fix the cutoff. 
Next the rest of the meson spectrum can be calculated using this cutoff: 

^ll^Lated = ^'{j\M'{A)\j), (9.2) 

where j refers to the state being calculated. 

^°In addition since the second C = — state is actually n — 4, we can not trust the accuracy of 
such an excited state. 

■^^We fix the lowest mass state for the given quantum numbers. 
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Secondly we analyze the spectrum for a range of a and ^ to find the values that 
give the correct values for the lowest and second lowest mass states above the state 
we fix. Finding which values of a and ^ correctly determine these two observables 
gives us our best values of a and ^ . 

9.3 Convergence Testing 

We must be certain we are using an adequate number of basis functions when 
taking final data. We verify our results have converged for various values of ^ with 
a = |.0 We reiterate that it is not necessary for the ground state eigenvalue to 
decrease as each additional function is added since the previous functions are shifted 

B 

Increasing the order of the B-splines (m in B-spline notation) beyond third order 
does not improve convergence, but increases the number of nonzero matrix elements. 
Therefore all B-splines we use are third order. Figures P?T] through pTPI show the 



convergence of the lowest three states for ^ = .2, .8, and 1.5 as we increase the number 
of longitudinal and transverse states. The number of B-splines in the longitudinal 
direction is: 

A;i + m + l-2 = A;i + 2, (9.3) 

where two states are not used because they have the incorrect behavior as x — 0, 1. 
We pair the splines into symmetric and antisymmetric functions under particle ex- 
change (see Sec. |6.4.1|) . Thus the total number of splines in the longitudinal direction 



^^We expect to need more longitudinal basis functions as ^ increases because the wavefunction 
should become more sharply peaked about '3^ — \ for heavy quarks. 

33c ~ 



^See Chapter H for details. 
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Figure 9.1: Convergence of the ground state for ^ = .2 and a = .5. 



Ni = ^ + 1. 
2 



Similarly the number of B-splines in the transverse direction is: 

Nt = ka + m + 1 - 3 = ka + 1, 



(9.4) 



(9.5) 



where three states are not used to ensure finite kinetic energy due to small x singu- 
larities. 

As increases a larger number of longitudinal states are needed for the eigenvalues 
to converge. Although we could have used fewer basis functions for smaller values of 



, all of our calculations use ki = 8 and /c2 = 5. Using equations (|9.4|), (p.5|), and 
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Figure 9.2: Convergence of the first exited state for ^ = .2 and a = .5. 
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Figure 9.3: Convergence of the second excited state for ^ = .2 and a = .5. 
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Figure 9.4: Convergence of the ground state for 
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Figure 9.5: Convergence of the first exited state for ^ = .8 and a = .5. 
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Figure 9.6: Convergence of the second excited state for 



.8 and a = .5. 
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Figure 9.7: Convergence of the ground state for ^ = 1.5 and a = .5. 
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n=1 ,m/A=1 .5 




Figure 9.8: Convergence of the first exited state for ^ = 1.5 and a = .5. 




Figure 9.9: Convergence of the second excited state for ^ = 1.5 and a = .5. 
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noting that we have four spin states, the number of basis states we use is: 

4x AT; xA^t = 4x5x6 = 120. (9.6) 

Since the Hamiltonian is Hermitian, at most 7260 of the 14400 matrix elements are 
unique. Since not all of the B-splines spatially overlap all of the other B-splines, there 
are 1680 elements that are zero. For our choice of ki and k2 we have reduced the 
number of elements that need to be calculated by 23%. 

9.4 Error Estimation 

There are three sources of errors in our meson spectra. Since we do not use an 
infinite number of basis functions we can not exactly determine the eigenvalues of our 
Hamiltonian. Based on our convergence testing, we can choose the smallest number 
of basis functions that produce converged eigenvalues to within two percent. These 
eigenvalues are of the Hamiltonian (IMO) , so a two percent error in the eigenvalues 
will produce approximately a one percent error in the masses. 

The second source of errors is statistical. Each matrix element is calculated by 
Monte-Carlo integration which suffers from statistical uncertainty. We calculate the 
error in the eigenvalues by allowing each matrix element to vary randomly within 
its gaussian distribution. This is done many times to generate a large number of 
eigenvalue sets that are used to calculate the average eigenvalue and the errors in the 
eigenvalues. The largest error allowed in any eigenvalue is two percent, giving a one 
percent error in the masses. 

The final and most significant source of error is in our approximation of the theory: 
limiting the approximate states to a quark-antiquark pair and truncating the Hamilto- 
nian at second order in the coupling. We cannot reliably estimate the error produced 
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by our approximation to the real state without first solving the meson spectrum with 
\qqg) states. We can estimate the error from approximating the Hamiltonian. 

Our Hamiltonian contains a contribution from the kinetic energy and second-order 
interactions, so we should expect errors of order g'^ (a^); however, for the heaviest 
mesons the errors in the masses should be very small because they are dominated by 
the kinetic energy. We can write 

+ =KE + V(l + e) (9.7) 

where 5 is the fractional error in the Hamiltonian and e is the fractional error in the 
interactions = g'^ (x a). However, what we actually solved is 

n = KE + Y. (9.8) 

We can now solve for the fractional error in the binding energy by considering the 
expectation values of Ti and KE: 

where we used the fact that the corrections to the interactions are 0{a). This means 
the error in the eigenvalues is roughly a times the fraction of the eigenvalue that 
comes from the interaction. 

9.5 hh Spectrum 

We calculate the hh spectrum to verify that our method gives similar results to 
0, ^, which uses nonrelativistic reduction (NR Reduction). The splittings in the bb 
system are so small that our method can not determine them in a reasonable amount 
of time. However the comparison is useful to make sure that we get similar results 
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State 


Experimental 


NR Reduction 


Fully Relativistic 


T 


9.460 


9.4 


9.64 ± .1 ± .1 ± .14 


Xbo 


9.860 


9.9 


N/A 34 


Xbi 


9.893 


9.9 


9.87± .1 ± .1 ± .14 


Xb2 


9.913 


9.9 


9.88±.1±.1±.14 




10.023 


N/A 


9.86±.1±.1±.14 



Table 9.1: Bottomonium masses in GeV from experimental results, nonrelativistic 
reduction and our fully relativistic calculation. The NR Reduction and Fully Rela- 
tivistic calculations use a = A and m = 4.9 GeV. The errors given for the relativistic 
calculation are from using a finite number of basis states, statistics and our second- 
order approximation to the Hamiltonian. 



for the quark mass, cutoff and meson masses. With the small splittings, our method 
also has little predictive power for this system. 

We use the values a = A and = 1.38 from and determine the cutoff to 



be 3.6 GeV by fixing the mass of Xbo (the lowest C = -|- state) to data. Table pTll 
shows our results (Fully Relativistic) along with the experimental results as well as 
those calculated using NR Reduction. 3.5% of mass of the Xbo is due to interactions. 
Thus we believe the errors in our bottomonium spectrum due to our approximate 
Hamiltonian to be about 1.4%. Our spectrum agrees with the other results within 
errors. This method passes the simplest test. 

The wavefunctions for the T^, Xbo, Xbi, and x^j, with a = A and j = are plotted 
in figures p.lO| , |9.11| , p.l2| , and |9.13| respectively. 
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This is the particle used to find the cutoff. 
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Figure 9.10: wavefunction. x is the longitudinal momentum fraction carried by 
one particle and k is the relative transverse momentum in units of the cutoff. 




l! 10 



Figure 9.11: Xfco wavefunction. x is the longitudinal momentum fraction carried by 
one particle and k is the relative transverse momentum in units of the cutoff. 
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Figure 9.12: Xh wavefunction. x is the longitudinal momentum fraction carried by 
one particle and k is the relative transverse momentum in units of the cutoff. 




Figure 9.13: Xb2 wavefunction. x is the longitudinal momentum fraction carried by 
one particle and k is the relative transverse momentum in units of the cutoff. 
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Cutoff{A) VS. a and m/A 
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Figure 9.14: The cutoff (MeV) versus tlie coupling and the charm quark; mass divided 
by the cutoff. 



9.6 cc Spectrum 

Our approximations may work; best for the cc system. We find a good represen- 
tation of the mass spectrum, but find unusual values for the coupling and cutoff. If 
we use the values of the coupling and ^ from we do not fit the spectrum within 
errors. 

We first generate eigenvalues for .1 < a < .9 and .1 < ^ < 1 for states with 
positive charge conjugation since only one of the five lightest states in the cc spectrum 
has negative charge conjugation. r]c is used to determine the cutoff which is plotted 



(in MeV) against a and ^ in Figure |9.14 



Using the corresponding value of A for a and ^ , we plot the mass of the first 



(n=l) and second (n=2) excited positive-charge-conjugation states in figures |9?T5| and 
9.16| respectively. 
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mass n=1 
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Figure 9.15: Mass of the first excited positive-charge-conjugation state, Xco as a 
function of a and ^ . The experimental value is 3417 MeV and each color band 
represents a mass range equal to two percent of the experimental mass. 



mass n=2 
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Figure 9.16: Mass of the second excited positive-charge-conjugation state, Xd as 
a function of a and ^ . The experimental value is 3510 MeV and each color band 
represents a mass range equal to two percent of the experimental mass. 
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Parameter Success Test 




Figure 9.17: This figure displays the success for choices of a and ^ at predicting 
the first and second excited state masses. We define a correct prediction to be within 
two percent of the actual value. The value is zero (purple) if neither mass is correctly 
determined, one (red) if one mass is correctly predicted and two (yellow) if both are 
correct. 



We find the values a and ^ that correctly produce the ratios of the first and 
second exited states by finding the region where the mass of the first and second 



excited states are predicted to within two percent. Figure |9.17| shows where both 
masses are predicted within the error. Unfortunately fixing these two masses is not 
sufficient to determine a and ^ because there is a range of a and ^ that correctly 
predict the spectrum. Also as is lowered the eigenvalues become smaller, increasing 
the cutoff. All of the combinations of a and ^ that successfully predict the spectrum 
produce a cutoff and/or coupling that is too large to be consistent with our 0{g^) 
approximation. If we consider a = | and ^ = .28, we get a cutoff of 3.8 GeV and 
the mass of the charm quark equal to 1.05 GeV. It is clear the errors generated by 
our approximations are too large to successfully predict the charmonium spectrum. 
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State 


Exp erimental 


NR Reduction 


Fully Relativistic 


Vc 


2.9798 


3.0 


N/A 35 




3.097 


3.0 


2.981 ±0.03 ±0.03 ±0.075 


Xco 


3.415 


3.5 


3.114 ±0.03 ±0.03 ±0.078 


Xci 


3.510 


3.5 


3.142 ±0.03 ±0.03 ±0.079 


Xc2 


3.556 


3.5 


3.145 ±0.03 ±0.03 ±0.079 



Table 9.2: Charmonium masses in GeV from experimental results, nonrelativistic 
reduction and fully relativistic calculations with a = .5 and ^ = .88. The errors 
given for the Fully Relativistic calculation are from using a finite number of basis 
states, statistics and our second-order approximation to the Hamiltonian. 



If we use a = .5 and ^ = .88 from |^ we find the spectrum in Table |9.2| and 
determine the cutoff to be 1.7 GeV, and the charm quark mass to be 1.5 GeV. 
5% of the rjc mass is due to interactions, making the error from the approximate 
Hamiltonian in our spectrum about 2.5%. The spin-averaged probability densities 



for the rjc, J/ip, Xco Xci with a = .5 and ^ = .88 are shown in figures pTT8| through 
9.21| , respectively. 



This was the particle used to find the cutoff. 
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Figure 9.18: r],. wavefunction. x is the longitudinal momentum fraction carried by 
one particle and k is the relative transverse momentum in units of the cutoff. 




Figure 9.19: J/ip wavefunction. x is the longitudinal momentum fraction carried by 
one particle and k is the relative transverse momentum in units of the cutoff. 
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Figure 9.20: Xcq wavefunction. x is the longitudinal momentum fraction carried by 
one particle and k is the relative transverse momentum in units of the cutoff. 




Figure 9.21: Xn wavefunction. x is the longitudinal momentum fraction carried by 
one particle and k is the relative transverse momentum in units of the cutoff. 
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CHAPTER 10 



Conclusions and Future 

The method we use to regulate and renormahze the hght-front QCD Hamiltonian 
does not remove any physics, add extra parameters or fundamentally limit itself 
to particular problems (like heavy quark systems). The renormalized, second-order, 
light-front Hamiltonian is reasonably successful at determining the glueball spectrum. 
This same approximation is much less successful at determining the heavy meson 
spectra. 

The success of the glueball calculation shows our method can be successful. It is 
clear that to more successfully determine (with predictive power) the heavy meson 
spectra we must either calculate the renormalized Hamiltonian to higher order in the 
strong coupling and/or allow for more complicated external states {qqg). In addition 
to these improvements using basis functions that are more sharply peaked should 
improve the accuracy of the heaviest meson system [bb) where the difference in the 
masses are fractions of a percent. 

We determined the light-front Hamiltonian to second order in the strong coupling 
while approximating the real states as a quark- antiquark pair. The results from the 
glueball and meson mass calculations show our method can reproduce experimental 
results with varying degrees of success. The renormalized Hamiltonian we determined 
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is the building block of higher order calculations and will be a necessary part of future 
calculations. 

It is unclear if B-splines are an efficient basis to approximate the real states. 
They reproduced the glueball spectrum with fewer basis functions than an earlier 
calculation but appear inefficient at representing the sharply peaked heavy meson 
wave functions. Although the full flexibility of the B-splines was not explored, it is 
likely more efficient to use basis functions better tailored to the problem at hand. 

The analytic calculation of the renormalized Hamiltonian requires a large time 
investment and great care. However the analytic result is not very useful without 
an effective and efficient method to numerically calculate the integrals the analytic 
calculation yields. We developed an algorithm to calculate these integrals without 
wasting time calculating any matrix element more accurately than needed. We also 
generalized the algorithm to run on a single processor or in a parallel processor en- 
vironment, including clusters of desktop machines. This is a significant development 
because the algorithm is independent of the details of the integrand so higher order 
analytic calculations can utilize it with very little change. 

We successfully added to the systematic determination of the renormahzed light- 
front Hamiltonian by including masses. We explored the effectiveness of B-spline basis 
functions and developed an algorithm for calculating the renormalized Hamiltonian 
matrix. 
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APPENDIX A 



Light-Front Conventions 



In this Appendix we list the hght-front conventions used in this dissertation. For 



a more general introduction see the light-front review by A. Harindranath |43 
A.l Coordinates 



X ■ y 
d ■ X 



^x^y 
Id+x- 



+ y^ — x\y\ 
+ \d^x^ — d]^x\ 





/x+ 


\ 




/O 2 











1 

2 





^\ 




x^ 






2 










1 











x^ = 


x^ 







-1 







2 







-1 







\ x'^ 


) 




^0 





-l) 












-1/ 



A. 2 Gamma Matrices 



1' 



r 




7" 



r 




r 



a I 
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A.3 Pauli Matrices 



^ =l,10j ^[z ) ^1,0-1 

A. 4 Projection Operators 



Ah 
A. 



2 

oo\_/io\ 

2\i J \ 2i ) ~ \ J 

I f -i\ f -2i \ _ f 
2\ i ojio oj~ioi 



Note: 



A+ + A_ = 1 A2 = A+ A2 = A_ 

A±A^ = 7+7+ = 7-7- = (7+)t = 7T 

^ A± 7°7+ = 7^7" 7^7+ = -7+7i 



A. 5 Lagrangian 



Cqcd = -^i^^.i^e'^'^ + m^)V^, (A.l) 



where 



Fr = a'^^r - d^'A^, - gA^^^A^r^^' (A.2) 



and 



T>^ij^{d^ + igAlTa)^l^. (A.3) 
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Greek indices are Lorentz indices, c's are color indices, a is the color index, repeated 
indices are summed over, and the /'s are the SU (iVc) structure constants. The Tq's 
are the SU (3) Gell-Mann matrices which obey the commutation relation [T^, T;,] = 

ifabcTc and are given by: 

/010\ /O -i 0\ /100\ /OOl 



1 



Ti 



1 

yo 

/Q -i 







fO-iO 




/I o^ 






i 


j T3 = i 


0-10 






i^O 




i^O Oy 






(0 0\ 




/O \ 


/ 




1 




-i 






lo 1 o) 




^0 i J 


V 



73 ^ ^ 

0^ 



A.6 The Gluon Field 

The expansion of the transverse gluon field takes the form: 
Ac±{x'',xl) = J DiSc,c, 



where 



V is the four-momentum operator, e is a positive infinitesimal, and 



The creation and annihilation operators follow the convention 
and have the commutation relations 



V3 / 



(A.4) 



(A.5) 



(A.6) 



(A.7) 



(A.8) 
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and 



[ai.aj] = [at, at] = 0, 



where 



The gluon field commutation relation is: 



x+=y+=0 

A. 7 Gluon Polarization Vector 



The gluon polarization vector is 



e+ = 0, e 



2q_L • e_L 



and 



so that 



E6'[(A)6r(A) = -^r, 

A 



where 7]^ = 7]^ = yf = Q and =2. Also 



We can write them out: 



^±1 



We also find: 

(7_L ■ = -(7_L ■ t-^'^ 



( \ 


\_ 

I 



V2 





1 



( \ 



1 

\-T2J 



a± ■ e_i = -cr_L • 
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A. 8 The Fermion Field 



The dynamical fermion degree of freedom is ipj^ and can be expanded in terms of 
plane wave creation and annihilation operators at = 0, 



(A.18) 



{x-,x±) = j Di6c,ci [bW+{k, oY^-'' + dxv\{k, o)e-'^-'\ , (A. 19) 



where the creation and annihilation operators are labeled the same as in Section |A.6 
The field operators satisfy 



x+=j/+=0 



and the creation and annihilation operators satisfy: 



A. 9 Dirac Spinors 



The dirac spinors a) and f (p, a) satisfy 



m) u{j), cr) = 0, + m) t>(p, cr) = 0, 



(A.20) 



(A.21) 



(A.22) 



and 



u{p,a)u{p,a') = -v{p,a)v{p,a') = 2m6ao 
u{p,a)'y'^u{p,a') = v{p,a)'j'^v{p,a') = 2p^6a 



(A.23) 
(A.24) 



u{p,a)u{p,a) = ;^ + m, ^ v{p,a)v{p,a) = ;^ — m. (A. 25) 
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We use the normalizations: 



(A.26) 



We can define: 



u+kx = Vk+ x'x v+kx = Vk+ x'x- 



(A.27) 



Use the eigenvalue equation: 



A+Xa = x'x 



(A.28) 



to find the eigenvector with eigenvalue 1. This leads to two solutions which we write: 

(A.29) 

We define xx as the upper 2 components of x'x- Finally, we discover that: 



Xi 

2 



/1\ 




Vo/ 



Xii 

2 



1 





«S2(P2) = -yJibl) ^^S2(P2) = \/piP2 Ssi,S2- 



(A.30) 
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APPENDIX B 



Details of the Combination of the Divergent Part of the 
Self-Energy and the Instantaneous Interaction Below the 

Cutoff 



The non-vanishing part of the instantaneous interaction below the cutoff is given 
32, 



'^e{x-y-t) 



S\,S2 ' 

g-A-*A|,^g-2A-4Aj.KA/K 



If we make the following definitions: 



S9 

Sl,S2 

F{x,y,q,p) = Biix)Biiy)Btiq)Btip)A}_^^^J(l)')A,.s,-sM 



xe 



we get: 



X^ J d^qd^pjl ^dx£ ^ dy ^ F{x,y,qJ). (B.2) 



(x - y) 



■^^This is an intermediate step in the calculation of the instantaneous interaction in Section 7.2.1 
that is not explicitly given. 
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Now make the change of variables: 



q-p Q = q + p 

^ Q+R ^ Q-R 

and find the Jacobian: 



This gives us: 

Q + R Q 



Defining: 



1 



= F(x,y, , ), 



we get: 



rx-t 1 

X2 / dy- -X3. 

Je [x - y) 

Integrating by parts, 

j udv — uv — j vdu 

with 



dy 

dv = V = — \og(x — y) u = 

x-y 
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gives: 



/x-e dX^ 
dy\og(x-y) — 



(B.9) 



Consider the first term: 



X2Xs\og{x - y)\y: 



-Xi / dx f d^Qd^N 

4 J2e J 



Q + ^¥^N Q-^x^yN 
xF[x, y, , ) log(a; -y)\y= 

-Xi / dx / d^Qd^R 

4 J2e J X — y 

, Q + R Q — R. . , I 
^F\x, y, — - — , — — ) log(a; - y) \y=^ 



Xi ' dx j d^qd^p 



2 ' 2 



x-y 

X F {x, y, q, p) \og{x -y)\y=, 
— Xi [ dx [ dqqdppd(t)d(j)' 

J2e J X — y 



xA*J(l>')Aj_,,_,MBi{x)Biiy)B,{q)M^^ 



^g-A-4A|.^g-2A-4Aj.KA,K 



log{x -y)\y= 



(B.IO) 



If we let 



X, = X,J dcpdcP'A*,_^^_JcP')A,_,^_,M 



s — 



P_ 



we get: 



eX^ I dx [ dqqdss— — Bi{x)Bv{y)Bt{q)Bti{s^/e) 

J2e J X — y 



^g-A-4A|.^g-2A-4Af.KA,K 



log{x -y)\y=,. (B.ll) 
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So now we need to consider this in the e — > hmit. 

m^{x — yY + es^(l — xY + — 2(1 — — y)^/esqcos{^) 



Afk = 



[1 - x){l - y){x - y) 



m?x'^ + 



— x) 

m?{x — yY + ts^x^ + e^q^ — 2xei sq cos{'^) 



xy{x - y) 



2 

2 



^FI 



m 

hs 

e 

+ es^ + q^ 



y{l - y) x{l - x) 
-^ + .^-^^. (B.12) 

The B-sphnes also go to zero as e — > 0. However, it is sufficient to note the exponential 
terms are the dominant factor in this limit, thus this term is zero. 
The second term is: 

X2X^\0g{x - y)\y=:c-e 

= /' ' dx I dqqdppd(t)d(t>'-\og{e)Bi{x)Bi,{y)Bt{q)Bt,{p) 

x^;'-..-..('/'')A-.i-..(0)e-^"^-e-^^-^^-^-|,=._,. (B.13) 

We can use the fact that when the angular integrals are done, we get 5^^/ and the 
change of variables; — — 0', to get: 



X2X3\og{x - y)\y=^ 

^"'-.X^ ['~'dx I dqqdppd(j)d<P'Bi{x)Biiy)Bt{q)Bt'{p) 

J2e J 



27r 

X (7,)e-^-^^-e-^^-^^-^- (B.14) 



e 

Next, when we integrate over 7 we get terms coming from the cosine in the angular 
basis function. We can also use the change of variables we used for the first term to 
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get: 



'-X^l^'^dx J dPQd'N\og{e)Bi{x)Bv{y)Bt{Q/2)Bt,{Q/2) 



Stt 

J 



X cos(L7- - - S2]70)e-^^^-e-^^-^^--^- (B.15) 
Use the following definitions: 

— * — * 

(5 = Q cos ax + (5 sin ay, N — N cosrjx + N sin-qy, (B.16) 



q = \\Jq^ + eiV2 + 2^eQNcos{a - rj), 
p = l\/Q^ + eN''- 2^/~eQNcos{a - rj), 



q-p- \{Q'-eN'). (B.17) 
With these, we know that: 

{p-^^^eN^ p'-q' = 0{^e), (B.18) 

Afk = -Aki ^ N^, Afi ^ 0. (B.19) 

Also, if we use: 

p-q — pqcos^, (B.20) 

we find that cos 7 = 1. This means that 7 = 0, so that cos{[j — si — 52)7) = 1. This 
leaves us with: 



^2^3log(a; - y)\y=x-e 

-^Xi ['~'dx [ d'Qd'Nlog{e)Bi{x)Bi,{x)Bt{Q/2)MQme~'''~"''" 
-^Xi / dx / d^Q\og{e)Bl{x)Bv{x)Bt{Q|2)Bt,{Q|2)t^ — ^f^. (B.21) 
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Next, if we do the angular integral in Q and let Q ^ 2q: 

J dFQf{Q/2) = 27r y dQQfiQ/2) ^SttJ dqqf{q), (B.22) 

we get: 



^2^3log(a: - y)\y=x-e 
= '"^^ / dqqlog{e)Bi{x)Bi,{x)Bt{q)BM 

= 5,yA2^(27r)t 1 cigglog(e)A(a;)A'(^)^t(?)^t'(?)- (B.23) 

So after adding the overall momentum-conserving delta function and the total longi- 
tudinal momenta which were dropped for convenience, we find: 

X,Xs \og(x - = -{q', I', t' ,j'\vi\A)\q, I, tj)"". (B.24) 

Thus, the divergence in the self energy is cancelled. 
The final term we need to calculate is: 



/x—e 
dy\og{x-y) 



dy 



dXs 
dy 



d^Qd^N dxj^" ' dy\og{x-y) 
= ^Xi ' dx r ' dylog(x-y)4- I d^Qd'^NX^ 

4 J2e Je dy J 

^Xi dx dy log{x -y)— / dV<^V 

J2€ Je dy J X — y 

rl—e rx—c r 

— Xi dx dy log{x — y) d^qd^p 

J2e Je J 



I2e 

xBi{x)BM)Mp)A],.,,^,M)^j-s.-sA<I>) 



x-^^S,,(y)e-^"'^-e-^^"'^--^-. (B.25) 
dyx — y 

We can rewrite A^/ in terms of A^^^ and A/^- After simplifying, we get: 
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Xi [ ' dx r ' dy\og{x-y) [ d^qd^p 

J2e Je J 



d 1 

dyx-y 



Bu{y)e 



-A-4(A|^+Af^) 



(B.26) 



or 



pl — c f-x—e r 

— Xi I dx dy log{x — y) d'^qd^p 

J2e Je J 



X 



J -S1-S2 



+ 



Kl) 



(B.27) 



{x-yY x-y x-y 

As we have seen before, when we do the sum over spins and do the integral over 0', 
we get: 

QO 



£ 'dxf^ ^dy\og{x-y) j d^qd^pBi{.x)Bt{q)Bt>{p)W^,^,e-''' 
" B,{y) , My) 



+ 



- 2- 



{x-yY x-y x-y 



Kl) 



, (B.28) 



where, 



Wg^q> = C0S(L7' - 1]7^) + 5g,2V,2COs([j + l]7<^) 

+ ((^9,3^,3 + (^9,4^,4) cos( j7^) , (B.29) 



m^(x — yY + (icp — yg)^ 

xy{x - y) 
m'^{x — y) — yq^ + xp^ {p — q) 



xy 



+ 



x-y 



^fk 



m^{x - yf + [(1 - x)p- (1 - y)qf 
{I - x){l - y){x - y) 
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m'^(x -y) + {l- y)q^ - (1 - x)p^ (p- q)' 



x){l-y) 



x-y 



A' 



KI 



{p — q)^ + J? 

(x - yY y2 



A' 



9 9 

m + p 



{p 



FK 



[l-yy {x-yf 

and we have left all six integrals intact. This means when we get to the end and 
are back in the original coordinates, we can just remove the integral over 0' as it has 
already been done. This integral appears divergent because if p = g, as x — > y the 
integrand diverges but the argument of the outer exponential is zero. The divergent 
part should cancel when we integrate over the transverse momenta. In the limit that 
X — y — T] : T] 0, the integrals over transverse momenta will bring down a factor 
of T]^. Thus the only term that may be divergent is the last one proportional to 
Afk^'fk + ^ki^'ki- However, we will consider the first and third terms: 
32 



-91y^JJ''P''H'P - dx dylogr) J d^qd^pBi{x)Bt{q)B,{p)W, 

^B,{y) My) 



- 2- 



T]^ r) 
Now consider writing this where x ~ y: 
32 n-e 



{Afk^FK + ^KlA'j^j) 



(B.30) 



-9lY^j,f'P''S('P-'P') dx dylogr) J d^qd^pBi{x)Bi,{y)Bt{q)Bt,{p) 



r]^ rj 

Note that in this limit, 7<^ 0, so that Wq^qi Sq^qi. Now let: 



:b.3i) 



1 q — p 
w — -- 



2 ' 



(B.32) 



where 



f—r [cos ax + sin ay] , w — w [cos 5x + sin t 
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(B.33) 



and 



(5 — a — 5, f-w — rw cos /3. 



(B.34) 



Then 



q — f+r)w, p — f—Tjw, 



(B.35) 



and 



q — ^Jr"^ -\- rjw'^ + 2rw^ cos /3, p = y^r^ + Tyty^ — Irw^ cos /3. (B.36) 



Define: 



r± = J r"^ + rjw^ ± 2rwy/rj cos /3. 



(B.37) 



Then 



2 2 

q ■ p — r — rjw , 



COS70 = , sin70 = — smp. (B.38) 



We have: 



(B.39) 



which gives: 



6 J2e Je V J 



1 - Uw^ 



(B.40) 



However, 



y dww (l - 64«;^) e-^^^'^ ^ q, 



(B.41) 



so the integral is finite. 
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Now let us try to rewrite the finite part of tlie apparent divergent integral in the 
original variables. 



1— e cx—t 



128 

= -gi—^3,j' dx I dyrjlogr] / d^rd^wBi{x)Bi,{y)Bt{r)Bt>{r)5q^q, 



1 - 64w^ 



32 /■! — e rx—e p _ „ 

-al^^j,]' J dyhgr] J d?qd^pBi{x)Bi'[y)Bt{r)Bt'{r)5y^q' 



32 



1— e px—e 



1-4 



-dAY^j'f I dyhgr] j d qd pBi{x)Bv{y)Bt{r)Bt'{r)6q^q. 



(B.42) 



Before the subtraction we had (Eq. |B.28|) 



(g',/',t',/|A^2(A)|g,/,t,j)f^ 



32 ("1 — e rx—t p 

-alY^jj'V^SiV - V) j^^ dx dy\og{x -y) j d^qd^p 



Bi'iy) , B'.iy) Biiy) 



+ 



- 2- 



(x — yy X — y X — y 



, (B.43) 



which when we change to the variables r and w gives: 



B,F 



{q'J',t',f\M\A)\qJ,t,j)^^ 



128 



1— e px—e 



-gl—Sjj'V^6{V -V) / dx / dyr]\ogr] / 



X 



+ 



- 2- 



?7 



(B.44) 



So after making the subtraction, we have: 
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{q',l',t',j'\M'iA)\q,l,t,j)^/^ 

128 r 
-gl^kf'P^SiV-V') I dxBiix) dyr]\ogr] j (frd 

My) 



w 



Bi'jy) B[,{y) 



rj 



Kl) 

1 - 64w^ 



?7 

-Bt{r)Bt,{r)Bi,{yy-^e 



77' 



Making the final changes of variables: 



and 



r] — xe 



-p 



dy = rjdp, 



(B.45) 



(B.46) 



1 + w 



1, 



1 + r 



1, 



1 +p 



1, 



(g',r,t',/|A<^(A)|g,/,t,j)fi^^ = 

pi p 27r /" 1 



where, 



dyr 



dy^ 



{l + y.fJ-i {l + yrYJ-i{l + y^f 
xr]logr) r wBi{x) 

+ B'M - 2B,{y) {AfkA'p^ + AkiA'^j) 



X 



(B.47) 



-B,{r)Mr)B,{y)^-^e-''-' 



1 - 64w' 



, (B.48) 



^q,q' = C0S([7' - 1]'^^) + 5<,,2(^g',2 COs([j + l]'y^) 

+ ((^9,3^,3 + <^9,4(^9',4) C0S( j7<^) , (B.49) 
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m^r) — yr\ + xr^ ^ ^ 2 



xy 



m 



'77 + (1 — y)r\ — (1 — x)r'i 



{l-x){l-y) 



rii'^ + ,.2 



V 



(B.50) 
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